5TH  SEMINAR 


BEST  ESTIMATE  OF  TRAJECTORY  TECHNIQUES 


2  OCTOBER  1978 


Space  &  Missile  Test  Center-Eastern  Test  Range 
Patrick  AFB,  Florida 


DATA  REDUCTION  AND  COMPUTING  GROUP 
RANGE  COMMANDERS  COUNCIL 


KWAJALEIN  MISSILE  RANGE 

WHITE  SANDS  MISSILE  RANGRmis  noaWI  IS  T»rsT  OmLTTY  PPACTICABI^. 

YUMA  PROVING  GROUND  THE  COPY  -  TO  TCC  COSTAIB®  A  / 

flBG'NIFI  CAJS X  HUJSER  OF  PAGES  WHICH  DO  -HO? 

NAVAL  WEAPONS  CENTER  HEFRODUCE  LEGIBLY., _ 

PACIFIC  MISSILE  TEST  CENTER 
ATLANTIC  FLEET  WEAPONS  TRAINING  FACILITY 
NAVAL  AIR  TEST  CENTER 


AIR  FORCE  FUGHT  TEST  CENTER 
AIR  FORCE  SATELLITE  CONTROL  FACILITY 
SPACE  AND  MISSILE  TEST  CENTER 
EASTERN  TEST.  RANGE 
WESTERN  TEST  RANGE 
ARMAMENT  DEVELOPMENT  AND  TEST  CENTER 
AIR  FORCE  TACTICAL  FIGHTER  WEAPONS  GENTBI 


Ttt  i  document  mM  bXB  arvnrovw!  o 
for  pubhc  rmloam*  and  oalr;  iis  ft 


dirfributtoa  Lb  unlimited. 


7«»li°08W4* 


DISCLAIMER  NOTICE 


THIS  DOCUMENT  IS  BEST  QUALITY 
PRACTICABLE.  THE  COPY  FURNISHED 
TO  DDC  CONTAINED  A  SIGNIFICANT 
NUMBER  OF  PAGES  WHICH  DO  NOT 
REPRODUCE  LEGIBLY. 


SECURITY  CLASSIFICATION  OF  THIS  PACE  (Whan  Data  Entmrmd) 


REPORT  DOCUMENTATION  PAGE 


1.  REPORT  NUMBER 


4.  TITLE  (end  Subtitle) 

5th  Seminar  -  Best  Estimate  of  Trajectory 
Techniques 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


3.  RECIPIENT’S  CATALOG  NUMBER 


5.  TYPE  OF  REPORT  ft  PERIOD  COVERED 


6.  PERFORMING  ORG.  REPORT  NUMBER 


Da1:a deduction  &  Computing  Group 

Range  Commanders  Council 

White  Sands  Missile  Range,  NM  88002 


9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Same  as  Block  7 


0.  contract  or  grant  numbers; 


10.  PROGRAM  ELEMENT.  PROJECT,  TASK 
AREA  ft  WORK  UNIT  NUMBERS 


It.  CONTROLLING  OFFICE  NAME  AND  ADDRESS  12.  REPORT  DATE 

Secretariat,  Range  Commanders  Council  2  October  1978 

ATTN:  STEWS-SA-R  is.  number  of  paces 

White  Sands  Msl  Rng,  NM  88002  280 _ 


4  MONITORING  AGENCY  NAME  ft  ADDRESS (if  different  from  Controlling  Office)  15.  SECURITY  CLASS,  (ot  thia  report) 


Same  as  Block  11 


UNCLASSIFIED 


15a.  OECLASSIFICATION/ DOWNGRADING 
SCHEDULE 


16.  DISTRIBUTION  STATEMENT  (of  thia  Report) 


This  docum>r:  hc~  1 
\m  public  -l  '  <x«.'5  ir- 

diatnl  >'.ti ••n  i.' 


17.  DISTRIBUTION  STATEMENT  (ol  the  abstract  entered  In  Block  20,  II  dltlerent  from  Report) 


19.  KEY  WORDS  (Continue  on  reverae  aide  if  neceaaery  end  identify  by  block  number) 

Best  Estimate  of  Trajectory  (BET),  black  boxes,  instrumentation/computer 
systems,  guidance  and  weapon  delivery  systems,  sophisticated  weaponry, 
accuracy  problems,  rigid  regression,  underspecified  models,  N  Interval 
Trajectory  Estimation  (NITE),  GPSBET ,  BET  development. 


20  ABSTRACT  (Continue  on  rereree  aide  If  neceaaery  end  Identity  by  block  number)  PREFACE 

\ie  Best  Estimate  of  Trajectory  (BET)  Seminar  was  conceived  when  the  Range 
Commanders  Council  Executive  Committee  and  the  Data  Reduction  and  Computing 
(roup  (DR&CG)  mutually  agreed  that  since  the  DR&CG  had  not  formally  addressed 
bET  in  several  years,  a  reexamination  of  this  technique  would  be  worthwhile. 

The  DR&CG  membership  felt  that  a  seminar  where  theoreticians  knowledgeable  in 
tii's  area  could  compare  and  contrast  BET  techniques  in  use  at  the  various 
nrges  would  provide  the  best  vehicle  for  bringing  new  material  on  this  subject 


EDITION  OF  I  NOV  «»  IS  OBSOLETE 


over 

UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 


UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  OF  THIS  PAOEQWnn  D«(«  Bnfrnd) 


to  the  surface. 


The  resulting  seminar,  which  was  held  on  2  October  1978,  was  a  valuable 
exercise.  All  speakers  agreed  that  although  theoretical  BET  problems  were 
solved  years  ago,  the  use  of  any  BET  program  still  requires  careful  analysis 
of  the  model  under  consideration.  Thus,  particular  attention  must  be  paid  to 
determining  if  enough  information  is  provided  to  solve  the  problem  at  hand  and 
whether  the  geometry  used,  combined  with  the  instrumentation  available,  will 
lead  to  an  ill-conditioned  problem.  It  was  noted  that  BET  programs  cannot 
presently  be  treated  as  "black  boxes."  However,  as  guidance  and  weapon  delivery 
systems  become  more  accurate,  the  search  for  instrumentation/computer  systems 
with  comparable  accuracies  that  can  be  used  as  a  standard  against  which  to  test 
this  sophisticated  weaponry  becomes  increasingly  more  challenging.  Since  BET 
techniques  comprise  in  most  cases  the  only  method  which  can  be  applied  to 
solving  sucfo  accuracy  problems,  their  importance  is  expected  to  increase 
appreciably  in  the  coming  years. 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  THIS  RAGEfW h*n  Dmm  t 


% 


TABLE  OF  CONTENTS 


PAGE 


LIST  OF  ATTENDEES .  v 

AGENDA .  ix 

PREFACE .  xi 

PAPERS 

NITE  -  N  Interval  Trajectory  Estimation  - 
Mr.  C.  W.  Welsh,  RCA,  Patrick  AFB,  FL .  1 

GPSBET  Program  Description  -  Mr.  R.  W.  Mai,  USAYPG  .  27 

TAPP  -  Trajectory  Analysis  and  Prediction  Program 
Overview  -  Mr.  G.  D.  Trimble,  Federal  Electric 
Corporation,  Vandenberg  AFB,  CA .  153 

BET  Development  at  WSMR  -  Mr.  W.  Agee,  WSMR .  203 

Ridge  Regression  with  Underspecified  Models  - 
Mr.  R.  Rowe,  RCA  Tech  Analysis,  Patrick  AFB,  FL .  243 


A  paper  entitled  "KMR  BET  Analysis,"  by  Tom  Keeney,  was 
not  received  in  time  for  publication. 


LIST  OF  ATTENDEES 


Name 

Danny  Weddle 

Robert  Barry 

William  A.  Thedford 

Richard  H.  Dale 

Wi 1 l iam  Agee 

William  L.  Johnson 

Richard  W.  Pace 

Jack  Harrison 

Edward  Gibson 


Address 


CS40/CSD 

Naval  Air  Test  Center 
Patuxent  River,  MD  20670 
AV:  356-3274 

6510  TESTW/TEESD 
AFFTC 

Edwards  AFB ,  CA  93523 
AV:  350-2240 

ANA-751 

DOT/FAA/NAFEC 

Atlantic  City,  NJ  08360 

CG 

STEWS-NR-AM 
ATTN:  R.  H.  Dale 
White  Sands  Missile  Range 
New  Mexico  88002 

CG 

STEWS-NR-AM 
ATTN:  W.  Agee 
White  Sands  Missile  Range 
New  Mexico  88002 

TFWC/PRDMX 

Nellis  AFB,  NV  89191 
CS42/CSD 

Naval  Air  Test  Center 
Patuxent  River,  MD  20670 
AV:  356-3365 

CSD 

Naval  Air  Test  Center 
Patuxent  River,  MD  20670 
AV:  356-4046 

ANA742 

DOT/FAA/NAFEC 
Atlantic  City,  NJ  08405 


r  ’ 

Name 

Address 

Shiu  M.  Cheung 

dot/faa/nafec 

ANA  751 

Atlantic  City,  NJ  08405 

Thomas  E.  Ford 

Code  3442 

Pacific  Missile  Test  Cente 
Point  Mugu,  CA  93042 

AV:  351-7931 

John  Greenwald 

Code  3153 

Pacific  Missile  Test  Cente 
Point  Mugu,  CA  93042 

AV:  351-6971 

1 1 

Paul  Welch 

ADTC 

Eglin  AFB ,  FL  32542 

AV:  872-3258 

B 

V.  B.  Kovac 

RCA-MTP 

T.  L.  Keeney 

BMDSC 

P.  0.  Box  1500 

Huntsville,  AL  35807 

Eugene  J.  Putzer 

Kentron,  Interra t i ona 1 

2003  Byrd  Spring  Road 
Huntsville,  AL  35802 

C.  W.  Welsh 

RCA 

MV  811,  Building  989 
Patrick  AFB,  FL  32925 

W.  M.  Nolin 

SAMTEC-ETR 

i 

G.  D.  Trimble 

PA300,  Building  6525 
Vandenberg  AFB,  CA  93454 
Ext.  7634 

i 

i 

D.  E.  Fletcher 

SAMTEC/ROCM 

Vandenberg  AFB,  CA  93437 
AV:  276-6746 

( 

! 

R.  N.  Learn 

NAVSWC 

Code  K-72 

Dahlgren ,  VA  22448 

AV:  249-7111 

VI 


Name 


Address 


J.  A.  Gilberto 


MAJ  Allen  Menard 


P.  Tokareff,  Jr. 


R .  R .  Rowe 


J.  A.  Greene 


Robert  W.  Mai 


Milton  Hillhouse 


Gerald  E.  Callahan 


J.  Robert  Gush 


USAOTEA 

CSTE-TDD 

5611  Columbia  Pike 
Falls  Church,  VA  22043 

TFWCRG/XRS 

Nellis  AFB ,  NV  89191 

RCA  Tech  Analysis 
DET-1  SAMTEC 
Patrick  AFB,  FL  32925 
MV  045 

AV:  494-5506 

RCA  Tech  Analysis 
DET-1  SAMTEC 
Patrick  AFB,  FL  32925 
MV  645 

AV:  494-4252 

RCA  Data  Processing 
DET-1  SAMTEC  PAFB ,  FL  32925 
AV:  494-7546 


Data  Reduction  Section 
STEYP-MDP 

Yuma  Proving  Grounds,  Yuma,  AZ  85364 
AV:  899-6767 

DEI  SAMTEC /TOOS 
Patrick  AFB,  FL  32925 
AV:  494-7735 

6510  Test  Wing 
AFFTC/TEEDOS,  STOPZOO 
Edwards  AFB,  CA  93523 

TFWCRG/XRSW 

Nellis  AFB,  NV  89191 


Vi  i 


AGENDA 


DATA  REDUCTION  AND  COMPUTING  GROUP 
OF  THE 
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BEST  ESTIMATE  OF  TRAJECTORY  (BET)  SEMINAR 
2  0CT02IR  1978 

MONDAY  -  2  OCTOBER  1978 

0800  OPENING  REMARKS  -  DANNY  WEDDLE,  NATC,  CHAIRMAN 

0805  NITE-N  INTERVAL  TRAJECTORY  ESTIMATION  -  C.  W.  WELSH, 

ETR/RCA 

0900  KMR  BET  ANALYSIS  -  TOM  KEENEY,  KMR 

1000  GPS  BET  -  ROBERT  W,  MAI,  YPG 

LUNCH 

1230  TAPP  -  TRAJECTORY  PREDICTION  PROGRAM  -  G.  TRIMBLE, 
SAMTEC/FEC 

1330  BET  DEVELOPMENT  AT  WSMR  -  BILL  AGEE,  WSMR 

1430  RIDGE  REGRESSION  WITH  UNDERSPECIFIED  MODELS  -  BOB  ROWE, 
ETR/RCA 

1530  CLOSING  REMARKS  -  DANNY  WEDDLE,  NATC,  CHAIRMAN 


PREFACE 


The  Best  Estimate  of  Trajectory  (BET)  Seminar  was  conceived  when 
the  Range  Commanders  Council  Executive  Committee  and  the  Data  Reduction 
and  Computing  Group  (DR&CG)  mutually  agreed  that  since  the  DR&CG  had 
not  formally  addressed  BET  in  several  years,  a  reexmaination  of  this 
technique  would  be  worthwhile.  The  DR&CG  membership  felt  that  a 
seminar  where  theoreticians  knowledgeable  in  this  area  could  compare 
and  contrast  BET  techniques  in  use  at  the  various  ranges  would  provide 
the  best  vehicle  for  bringing  new  material  on  this  subject  to  the 
surface. 

r 

The'resulting  seminar,  which  was  held  on  2  October  1978,  was  a 
valuable  exercise.  All  speakers  agreed  that  although  theoretical  BET 
problems  were  solved  years  ago,  the  use  of  any  BET  program  still  re¬ 
quires  careful  analysis  of  the  model  under  consideration.  Thus, 
particular  attention  must  be  paid  to  determining  if  enough  information 
is  provided  to  solve  the  problem  at  hand  and  whether  the  geometry  used, 
combined  with  the  instrumentation  available,  will  lead  to  an  ill- 
conditioned  problem.  It  was  noted  that  BET  programs  cannot  presently 
be  treated  as  "black  boxes."  However,  as  guidance  and  weapon  delivery 
systems  become  more  accurate,  the  search  for  instrumentation/computer 
systems  with  comparable  accuracies  that  can  be  used  as  a  standard 
against  which  to  test  this  sophisticated  weaponry  becomes  increasingly 
more  challenging.  Since  BET  techniques  comprise  in  most  cases  the  only 
method  which  can  be  applied  to  solving  such  accuracy  problems,  their 
importance  is  expected  to  increase  appreciably  in  the  coming  years. 
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NITE  -  N  INTERVAL  TRAJECTORY  ESTIMATION 
C.  W.  Welsh 

RCA,  Patrick  AFB ,  Florida 


ABSTRACT 


NITE  is  the  "Best  Estimate  of  Trajectory"  computer  program  developed 
at  SAMTEC  Det  #1  (ETR) ,  Patrick  AFB,  Florida.  NITE  is  a  minimum  variance 
(weighted  least  squares)  trajectory  and  error  model  coefficient  estimation 
program  designed  to  process  tracker  and  guidance  data  from  single  or  multi¬ 
ple  powered  flight  and/or  free  fall  trajectories  from  launch  to  impact  in 
a  single  computer  run.  Tracking  instrumentation  may  be  either  earth-fixed 
or  referenced  to  moving  ships,  aircraft  and  satellites  whose  trajectories 
may  also  be  simultaneously  estimated  (differentially  corrected).  The  pro¬ 
gram  will  accept  either  preprocessed  or  raw  data  and  will  perform  a  variance 
Geodetic  Dilution  of  Precision  (GDOP)  analysis  either  as  a  primary  task  or 

as  a  byproduct  of  the  trajectory  and  error  model  coefficient  estimation. 

The  variance  propagation  considers  both  estimated  (modeled)  and  enforced 
(unmodeled)  systematic  error  terms. 

NITE  has  been  used  to  provide  trajectory,  instrumentation  error 
model  coefficient,  station  location,  ship  location,  aircraft  location, 
drag  and  lift  coefficients  for  a  broad  class  of  launch  vehicles,  reentry 
vehicles  and  satellites  such  as  Minuteman  III,  Titan  ill,  Polaris, 

Poseidon,  Trident,  Thor-Delta,  Atlas-Centaur,  Saturn,  Pershing,  SRAM, 

Geos  A,  Geos  B,  Pageos,  Transit  and  Echo.  It  is  presently  being  used 
to  process  data  at  SAMTEC -Western  Test  Range  (WTR). 

NITE  was  originally  written  for  the  IBM  7094  computer  in  Fortran  IV 
with  a  special  IBSYS  monitor.  It  lias  also  been  modified  for  operation 
on  the  IBM  360/65  using  a  modified  OS  monitor,  and  for  operation  on  the 
CDC  Cyber  74  aixl  GDC  6500  using  the  NOS  monitor. 
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I .  INTRODUCTION 


NITE  is  an  acronym  used  for  the  N-lnterval  Trajectory  Estimation 
computer  program  to  produce  the  "Best  Estimate  of  Trajectory"  (BET)  at 
the  Eastern  Test  Range.  The  N-Intcrval  concept  evolved  from  the  desire 
to  combine  several  missile  and/or  satellite  trajectories  into  one  solu 
tion  to  enhance  the  estimation  capability  of  surveys  and/or  tracker 
error  model  coefficients.  However,  most  of  the  BEl"s  produced  at  El’R 
are  of  the  single  interval  type.  Sufficient  flexibility  has  been  designee', 
into  the  program  to  allow  the  use  of  earth-fixed  and/or  moving  trackers, 
inertial  and/or  earth  fixed  data,  with  a  wide  variety  of  types  of  measure¬ 
ments.  The  locations  of  all  (moving  or  fixed)  crackers  may  be  simulta¬ 
neously  estimated  (differentially  corrected).  The  input  measurements  may 
be  either  preprocessed  or  raw  (decoded)  data.  In  the  latter  case,  NITE 
will  correct  data  for  refraction,  transit  time,  scaling  and  a  compre¬ 
hensive  error  model.  NITE  will  also  perform  a  Variance  or  Geodetic  Dilute  n 
of  Precision  (GDOP)  analysis  as  a  primary  product  or  as  a  byproduct  of 
the  trajectory  arid  error  model  coefficient  estimation.  The  variance 
propagation  considers  both  the  estimated  (modeled)  and  enforced  (un;.  -.le 
systematic  error  terms. 

NITE's  prime  task  has  been  to  provide  trajectory,  instrumentation 
error  model  coefficient,  station  location,  ship  location,  aircraft  local; 
lift  and  drag  coefficient  estimates  for  all  major  ETR  launches.  NITE  i ; 
also  used  at  WTR. 

The  mathematical  basis  for  NITE  is  minimum  variance  (weighted  least 
squares)  which  can  be  found  in  most  statistics  textbooks.  References 
(2)  through  (8)  contain  related  discussions  of  trajectory  and  error  model 
coefficient  estimation  theory.  The  program  document  (1)  gives  a  detailed 
description  of  the  mechanics  of  the  KITE  estimation  and  error  analysis 
techniques.  Only  enough  mathematics  necessary  for  describing  instrumen¬ 
tation  error  models  and  for  interpreting  results  will  be  presented  in 
this  paper. 


The  design  of  NITE  is  the  result  of  many  interacting  forces.  lYimary 
among  these  forces  has  been  the  experience  gained  by  ETO  personnel  in 
processing  and  analyzing  data  from  thousands  of  missile  launches  aixi 
satellite  passes  over  the  past  18  years  and  in  the  development  of  earlier 
BET  programs  such  as  BEl'A,  BELY,  GLAD,  TRkD  and  TEEM.  Another  significant 
influence  in  the  program's  design  has  been  the  many  requirements  and 
suggestions  from  RCA  Technical  Analysis,  PAA  Tech  Staff,  TRW  Systems, 

The  Aerospace  Corporation,  NASA  and  from  RCA  and  PAA  consultants.  Many 
similar  computer  programs  developed  by  other  agencies  have  been  studied 
for  worthwhile  ideas  applicable  to  ETR  requirements.  The  programs  studied 
include  those  developed  by  the  PAA  Tech  Staff,  The  Aerospace  Corporation, 
JPL,  Lockheed,  TOW  Systems,  The  Mitre  Corporation,  and  The  Wolf  Research 
and  Development  Corporation  bearing  the  titles  ERAN,  ORAN,  TRACE,  REMP, 
TRAP,  EMBET,  DPODP  and  MARLOCK. 

NITE  was  originally  written  for  the  IBM  7094  computer  in  Fortran  IV 
using  a  special  IBSYS  monitor.  It  was  later  modified  for  operation  on 
the  IBM  360/65  using  a  special  OS  monitor.  It  has  been  recently  modified 
to  run  on  the  GDC  6500  and  CDC  CYBER-74  using  the  CDC  NOS  monitor. 
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11.  THEPURPOSG  01  THU  NITE  PROGRAM 


Lhe  original  purpose  of  the  EfR  BE!  programs  was  to  combine  all  of 
the  tracker  data  from  a  missile  launch  to  provide  the  range  user  with 
better  trajectory  estimates  than  those  from  a  single  tracker.  Although 
this  is  still  the  primary  goal ,  the  BE!  concept  has  been  extended  to 
cover  many  other  tasks  such  as  the  following: 

A.  fo  estimate  (differentially  correct)  multiple  powered  flight,  orbit. ll . 
sub-orbital  and  reentry  trajectories  in  a  single  run  while  simulta 
neously  estimating 

1.  coefficients  of  large  tracker  and  inertial  guidance  instrumen¬ 
tation  error  models; 

2.  geodetic  locations  of  eartn  fixed  trackers; 

3.  trajectories  of  shipboard,  airben.  ,  and  satellite-borne  tna-  r 

4.  drag  and/or  lift  error  model  coefficients; 

5.  geopotential  coefficients*, 

6.  impact  points; 

7.  observation  (measurement)  weights, 

8.  underwater  transponder  locations; 

/  9.  re  la  ive  ]  itions  ind  ve  ,  t,  f  multiple  reentry  bo 

10.  other  quai  vn<  se  partial  derivatives  can  be  computed  ext< 

nally  such  as 

a.  missile  pitch,  roll,  yaw  and  thrust  decay  coefficients. 

b.  ship  or  aircraft  pitch,  roll,  flexure  and  inertial  naviga; i>  n 
coefficients; 

c.  atmospheric  pressure,  electron  density,  tempt* rature ,  wind 
velocity; 
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d.  point  mass  geopotential  terms; 

e.  second  order  instrumentation  error  model  terms. 

This  capability  often  eliminates  the  need  for  developing  a  special 
analysis-type  program  when  a  new  tracker  system  or  a  different  physical 
phenomenon  requires  investigation. 

B.  To  perform  a  variance  (GDOP)  analysis  for  the  quantities  being  estimated 
or  unmodeled,  either  as  a  byproduct  of  the  estimation  process  or  as 

a  product  of  a  special  GDOP  run, using  theoretical  trajectories  and 
the  theoretical  trackers.  The  latter  type  runs  are  of  value  in  range 
planning. 

C.  To  provide  intersystem  data  comparisons. 

D.  To  perform  orbital  prediction  accuracy  analysis. 

E.  To  show  the  effect  of  ignoring  selected  error  model  terms  during  the 
estimation  process  (using  unmodeled  error  analysis). 

F.  To  apply  various  constraints  on  or  between  quantities  being  estimated. 

G.  To  correct  measurements  for  known  errors  such  as  refraction,  transit 
time  and  erroneous  calibrations. 
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HI.  NITF  PROGRAM  CAPABI LITIES 


Measurement  Types 

L.  Azimuth  position,  velocity  and  acceleration. 

2.  Elevation  position,  velocity  and  acceleration. 

3.  X-Y  mount  position,  velocity  and  acceleration. 

A.  Range  position,  velocity  aixi  acceleration. 

3.  Range  sum  position,  velocity  and  acceleration. 

6.  Range  difference  position,  velocity  and  acceleration. 

7.  Inertial  guidance  position  and  velocity. 

8.  Ballistic  camera  plate  coordinates. 

9.  iopocentric  and  geocentric  right  ascension  and  declination. 

10.  Earth  fixed  cartesian  coordinate  position  and  velocity  components. 

11.  Latitude,  longitude  and  height. 

Modes  of  Operation 

1 .  EDIT  Mode 

Input  observations  may  be  compared  with  a  transformed  reference 
trajectory  and  may  be  corrected  for  most  types  of  systematic 
errors.  The  corrected  output  of  the  EDIT  mode  run  may  be  used 
directly  as  input  to  an  estimation  mode  run. 

2 .  ESTIMATION  Mode 

A  minimum  variance  (weighted  least  squares)  estimate  of  the 
trajectory  and  error  model  coefficients  is  computed  from  the 
input  observations.  GDOPs,  residuals  and  partial  derivatives 
are  byproducts  of  the  estimation  mode. 


Coefficients  of  instrumentation,  drag,  lift,  geopotential ,  or 
other  error  models  may  be  estimated  or  may  be  unmodeled.  Un¬ 
modeled  is  a  term  used  for  the  process  that  applies  known  co¬ 
efficients  (not  estimated)  and  propagates  a  given  uncertainty. 

3.  SIMULATION  Mode 

Random  and  systematic  errors  may  be  applied  to  measurements  computed 
from  a  reference  trajectory.  The  output  of  the  simulation  mode  may 
be  used  as  input  to  estimation  or  edit  mode  runs.  Tnus  the  ability 
of  the  program  to  estimate  known  systematic  errors  under  various 
conditions  may  be  estimated. 

4.  ORBITAL  PREDICTION  Mode 

Free-fall  initial  conditions  are  determined  from  data  during  early 
satellite  passes.  The  trajectory  is  then  generated  and  compared 
with  data  acquired  on  the  same  satellite  many  passes  later. 

5.  GDOP  Mode 

Only  the  error  propagation  features  of  the  program  are  utilized. 
Estimates  of  variance  of  the  input  observations  (including  a  priori 

estimates  of  variance  in  survey  and  other  error  model 
terms  as  well  as  in  instrumentation  observations)  are  transformed 
into  estimates  of  variance  in  the  parameters  being  estimated.  The 
magnitude  of  the  GDOPs  is  solely  a  function  of  the  assumed  tra¬ 
jectory,  the  assumed  instrumentation  combination,  assumed  random 
observation  errors  and  the  a  priori  knowledge  of  the  systematic 
error  components.  The  GDOP  mode  may  thus  be  used  to  evaluate 
hypothetical  instrumentation  systems  operating  under  hypothetical 
conditions. 

Data  Handling 

1.  Data  may  be  on  either  tape  or  cards  or  both. 

2.  Data  may  be  on  one  or  several  files.  Tapes  will  be  automatical ly 
rewound  if  necessary. 
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k.  Dynamic  lag 

l.  Radar  mislevel ,  col  1 imat i on,  nonorthogonality,  drcxip 
in.  Tracker  survey  error  terms  of  severaL  types 

(1)  Movement  of  latitude,  longitude  and  weight 

(2)  Movement  ol  X,  Y,  L  rectangular  cartesian  coordinates 

(•1)  1 rans 1 at ion  only  with  fixed  rotation  (for  astronomically 
ori ent  ed  systems ) 

(4)  Both  translation  and  rotation 

(5)  Movement  of  a  central  site,  with  related  stations  moving 
as  an  array 

(6)  Base  line  azimuth,  elevation  and  range  from  the  central  site 

(7)  Combinations  of  items  1-6,  allowing  both  "external"  and 
"internal"  survey  adjustments 

n.  Tracker  velocity  if  on  moving  aircraft,  ship  or  satellite 

(1)  latitude,  longitude  and  altitude  rates 

(2)  Heading  and  speed 

o.  Random  errors  for  simulation 

3.  Products  of  terms  may  be  formed.  This  allows  the  following 
additional  partials: 

a.  The  product  of  a  timing  error  partial  and  a  reinitialization 
partial,  allowing  correction  for  a  timing  error  that  occurs 
within  a  data  span. 

b.  The  products  of  survey  partials  and  time  polynomial  partials 
allow  ship  or  aircraft  movement  to  be  expressed  as  poly¬ 
nomials  in  time. 

c.  The  product  of  trigonometric  functions  aLlowing  modeling  of 
second  order  radar  errors. 

d-  The  products  of  values  from  the  input  tape  and  built-in  partials 
can  be  used  to  expand  instrumentation  error  models. 
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4.  The  need  for  estimating  a  zero  set  term  in  free  fall  can  be 
eliminated  while  still  considering  its  effect  on  the  solution 
by  defining  a  measurement  as  a  reinitialization  type.  The 
measurement  is  then  defined  as  a  difference  between  the  current 
value  and  the  value  at  the  first  time  the  measurement  entered 
the  solution.  This  feature  was  designed  mainly  for  long  arc 
work,  but  can  also  be  used  for  SRN-9  system  ship  location. 

Free-Fall  Long  and  Short  Arc  Capabilities 

1.  Up  to  n  short  arc  initial  conditions  may  be  estimated  (where  n 
is  the  number  of  intervals:  unlimited  for  the  7094  version  of 
NITE,  limited  to  14  for  the  360/65  version  and  limited  to  9  for 
the  CYBER/ 6500  version).  A  long  arc  (a  long  arc  is  defined  as 
one  that  extends  over  more  than  one  interval;  see  F-2  below) 
may  be  preceded  by  any  number  of  short  arcs  or  powered  flight 
trajectories  but  must  be  the  last  arc  estimated  in  the  run. 

2.  Fen  different  body  types  defined  by  weight,  cross-sectional  re¬ 
flection  area,  coefficients  of  drag,  lift  and  reflectivity  may 
be  handled  in  one  computer  run.  A  drag  table  (mach  number  vs 
drag  coefficient,  altitude  vs  drag  coefficient  or  time  vs  drag 
coefficient)  may  be  applied  to  each  object.  Lift  tables  may 

be  entered  in  the  same  manner.  Up  to  300  entries  are  allowed 
in  each  table  in  the  360/65  and  CYBER/6500  versions  of  NITE. 

3.  Any  desired  geopotential  model  may  be  read  in  on  cards.  The 
maximum  size  is  order  24  degree  24  for  the  7094  and  360/65 
version  and  order  30  degree  30  for  the  CYBER/6500  version. 

ihe  7094  version  has  6  built-in  geopotential  models,  the  360/65 
version  has  21  built-in  geopotential  models  and  the  CYBER/6500 
version  has  5  built-in  geopotential  models. 

A1 1  models  may  be  truncated  by  the  program  user.  Up  to  15  terms 
may  be  estimated  or  modified  through  the  station  constants,  ihe 
built-in  models  are  updated  as  improved  models  become  available. 
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4.  Up  Lo  15  of  the  constants  associated  with  free  fall  (drag,  lift, 
cross-sectional  area,  geopotential  coefficients,  etc.)  may  be 
estimated  or  have  their  a  priori  variances  propagated  into  the 
GDQPs  on  one  computer  run. 

5.  The  Cowell  Method,  using  a  Gauss- Jackson  with  automatic  step 
size  selection,  is  used  Ln  orbital  integration,  Ihe  integration 
is  eighth  order. 

6.  ihe  built-in  atmospheric  model  is  the  Patrick  reference  below 
30  km.,  U.S.  Standard  (1962)  Model  between  30  km  and  120  km, 
and  a  modified  Jacchia  1964  (Dynamic)  model  above  120  km. 
Atmosphere  models  may  also  be  read  in  through  the  station  con¬ 
stants  either  as  card- by- card  samples  (standard  weather  deck) 
or  as  polynomials.  Wind  speed  and  wind  direction  profiles  may 
also  be  entered. 

7.  Extraterrestrial  potential  (lunar,  solar  and  planetary)  and 
solar  radiation  pressure  may  be  applied  to  all  objects. 

8.  Initial  conditions  (initial  position  and  velocity  vector)  used 
as  input  may  be  of  10  different  types. 

9.  A  self- start  capability  allows  initial  conditions  to  be  computed 
from  two  position  vectors. 

10.  Six  types  of  orbital  elements  are  output  options. 

11.  Discontinuities  in  position  and  velocity  (velocity  kicks)  may 
be  entered  through  the  error  model  coefficients. 

12.  Trajectories  may  be  constrained  to  pass  tnrough  known  impact 
points  or  to  impact  at  specified  times. 

Number  of  Parameters  and  Measurements  A1,  lowed 

1.  By  building  up  the  normal  equation  coefficients  one  at  a  time  and 
by  making  use  of  the  fact  that  element  i,  i  of  the  normal  equation 
coefficient  matrix  is  zero  if  parameter  i  is  never  in  the  solution 
at  the  same  time  or  during  the  same  arc  as  parameter  i,  the  only 
limit  to  the  number  of  measurements  instrumentation  error  model 
terms,  trajectory  coordinates  or  free-fall  initial  conditions 
which  may  be  estimated  during  a  computer  run  is  the  number  that 
can  be  defined  in  the  n  intervals. 
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3.  Data  from  one  NITE  run  may  be  used  as  input  to  other  NITE  runs. 

4.  Data  may  be  eliminated  from  solutions  through  start -stop  time 
control  by  elevation  angle  cutoffs,  or  through  a  missing  data 
indicator  (all  bits  for  the  input  data  tape  word). 

5.  Data  may  be  generated  from  either  the  input  or  output  trajectories. 
Error  Model  Terms 

1.  Any  terms  whose  partial  derivatives  can  be  computed  externally 
and  placed  on  the  NI'IE  input  tape.  Inertial  guidance  terms, 
many  of  which  require  numerical  integration, have  been  handled 
in  this  manner.  Complex  atmospheric  models,  thrust  models  and 
second-order  instrumentation  error  models  may  also  be  handled 
in  this  way. 

2.  The  following  terms  are  built  into  the  NITE  partial  derivative 
subroutine : 

a.  Zero  sets  for  all  the  instrumentations 

b.  Reinitialization  of  measurements  at  specified  times.  This 
allows  compensation  for  shifts  in  the  input  measurements  in 
addition  to  the  original  bias. 

c.  Velocity  measurement  biases 

d.  Acceleration  measurement  biases 

e.  Rate  biases  caused  by  an  error  in  the  local  frequency  reference 

f.  Polynomials  in  time  for  trajectories  and  all  input  measurements 

g-  Scale  factor  (to  compensate  for  errors  in  wave  length,  speed 
of  light  or  sound,  focal  length,  etc.) 

h.  Refraction  correction  of  the  range,  range  sum,  range  difference, 
elevation  and  ballistic  camera  plate  coordinates 

i  •  Timing  biases 

j .  Transit  time  corrections 
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2.  In  order  to  make  efficient  use  of  the  sparse  nature  of  the  error 
model  normal  equation  coefficient  matrix,  the  program  breaks  a 
computer  run  into  intervals  during  which  specified  error  nxxiel 
terns  do  not  occur  simultaneously  in  the  solution  with  other 
terms.  This  allows  the  inversion,  by  partitioning,  of  a  matrix 
of  an  order  limited  only  by  the  number  of  intervals.  The  largest 
number  of  error  model  terms  which  can  be  handled  during  any 
interval  is  100.  If  m  specified  terms  (such  as  survey  and  geo¬ 
potential)  are  in  the  solution  during  more  than  one  interval, 
100-m  additional  terms  may  be  estimated  in  eacn  of  the  intervals. 
The  maximum  number  of  terns  which  can  be  estimated  during  a 
computer  run  broken  into  n  intervals  would  then  be  m+n  (100-m) . 
Thus,  if  each  term  were  in  the  solution  during  more  than  one 
interval,  the  limit  for  the  computer  run  would  be  100  terms. 

3.  Three,  six  or  nine  (position,  position  and  velocity,  position 
and  velocity,  and  acceleration)  trajectory  parameters  may  be 
estimated  at  each  timepoint. 

4.  Forty  position,  40  velocity  and  40  acceleration  type 
measurements  may  be  used  at  each  time,  and  in  each  interval . 

Thus,  if  there  were  m  times  per  interval,  n  intervals  and  40 
position,  40  velocity  and  40  acceleration  measurements 

at  each  time,  the  total  number  of  measurements  for  the  run 
would  be  120xnxm. 

Other  Features 

1.  The  weighting  of  measurements  may  be  performed  in  several  wavs: 

a  .  Weights  may  be  input  from  tape 

b.  Weights  may  be  input  by  station  constants 

c  .  Weights  may  be  computed  by  the  program  as  functions  of  error 
model  terms 

d.  Any  combination  of  a,  b.  and  c.  may  be  used 
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2.  A  full  a  priori  covariance  matrix  for  the  first  ten  instrumentation 
error  models  may  be  input. 

3.  Long  running  7094  computer  runs  may  be  interrupted  and  restarted 
at  a  later  time.  This  is  accomplished  by  saving  all  information 
accumulated  up  to  the  time  of  interruption.  The  360/65  version 
has  been  adapted  to  run  from  a  CRT  terminal ,  where  data  are 
displayed  on  the  CRT  during  the  run.  The  CYBER/6500  version 
can  also  be  run  from  a  terminal . 

4.  Dual  tracking  points  may  be  assumed  for  range  sum  and  range 
difference  type  measurements. 

5.  Tracking  points  may  be  switched  at  specified  times.  This  is 
sometimes  necessary  when  different  antennas  are  used  after 
staging  events. 

6.  Attitude  information  necessary  in  computing  tracking  point  cor¬ 
rections  may  either  be  derived  from  velocity  vectors  or  input 
as  pitch,  roll  and  yaw  from  tape.  Attitude  rate  data  used  in 
correcting  velocity  data  may  also  be  input. 

7.  Errors  in  parameter  estimates  due  to  unit  errors  in  unmodeled 
error  coefficients  may  be  computed  as  an  option. 

8.  Improved  estimates  of  instrumentation  measurement  weights  may 
be  computed  from  residuals  as  an  option. 

9.  Both  convergence  and  divergence  tolerances  are  tested  in  order 
to  prevent  unnecessary  calculations. 

10.  Partials  of  residuals  with  respect  to  unmodeled  errors  may  be 
output . 

11.  Residual  analysis  output  includes  means,  standard  deviations 
and  ranges  of  residuals. 

12.  Reentry  data  may  be  output:  altitude,  mach  number,  atmospheric 
density,  coefficients  of  drag  and  lift,  lift  vector  angle, 
ballistic  coefficient,  drag  and  lift  accelerations,  and  reentry 
angle. 
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13.  Ship  pitch,  roil  and  heading  and  their  velocities  may  be  input 
to  correct  shipboard  or  airborne  position  and  velocity  measure¬ 
ments. 

14.  Ships  oi-  aircraft  may  be  constrained  to  follow  either  great 
circle  or  rhunb  line  (constant  heading). 

15.  Plots  of  the  residuals  may  be  output  on  printer. 
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IV.  INPUT  REQUIREMENTS 

The  input  requirement  can  be  explained  by  following  the  structure  of  the 
station  constants,  which  are  segmented  into  groups. 

A.  Group  1  -  General  constants  common  to  all  intervals 

These  are  constants  for  decision  logic  such  as  the  number  of  intervals; 
outer  iteration  divergence  tolerance;  position  only,  position  and 
velocity,  or  position  and  velocity  and  acceleration;  mode  of  run  (GDOP, 
EDIT,  SIMULATION  or  ESTIMATION);  output  in  feet  or  meters,  right 
handed  or  left  handed;  maximum  number  of  words  to  read  from  the  input 
tape;  outer  iteration  convergence  tolerance;  number  of  outer  iterations; 
moving  origin  differencing;  frequency  of  printed  output;  frequency  of 
input;  number  of  inner  iterations;  and  EDIT  tolerances. 

B.  Group  2  -  Master  Coordinate  System 

This  defines  the  location  of  the  output  trajectory,  its  orientation 
and  spheroid. 

C.  Group  3A  -  Independent  and  Adjustable  Auxiliary  Origins 

Up  to  10  surveys  can  be  defined  in  this  group  by  geodetic  and 
astronomic  coordinates. 

D.  Group  3B  -  Survey  Covariance  Matrices 

Survey  covariance  matrices  may  be  defined  for  the  surveys  in  Group  3A. 

E.  Group  4A  -  General  Free-Fall  Constants 

This  group  defines  constants  (logical  and  data)  for  the  operation  of 
the  free-fall  constraints. 

F.  Group  4B  -  Estimation  of  Free-Fall  Parameters 

These  are  the  constants  to  define  up  to  15  free-fall  parameters  to 
be  estimated,  their  a  priori  value  and  their  a  priori  standard  deviation. 
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G.  Group  5A  -  Adjustable  System  Parameters 

Similar  to  4B,  up  to  100  adjustable  system  parameters  are  defined, 
and  their  a  priori  values  and  standard  deviations  are  given.  Any 
systematic  parameters  to  be  used  for  weights  are  defined,  as  are 
those  parameters  that  will  not  be  estimated,  but  whose  uncertainties 
will  be  propagated. 

H.  Group  5B  -  Covariance  of  Parameters 

The  covariances  of  up  to  10  parameters  defined  in  Group  5A  may  be 
input . 

I.  Group  6A  -  Independent,  Nonadjust able  Auxiliary  Origins 

Up  to  20  geodetic  and  astronomic  origins  may  be  defined  for  all  types 
oi  systems  including  inertial  systems.  Aspect  angles  relative  to 
particular  surveys  may  be  called. 

J.  Group  6B  -  Dependent,  Nonadjustable  Auxiliary  Origins 

Up  to  20  origins  may  be  defined  that  are  related  to  Group  3A  or 
Group  6A  origins  by  a  fixed  cartesian  coordinate  relationship. 

K.  Group  6C  -  Camera  Orientation 

Up  to  20  ballistic  camera  orientations  may  be  defined. 

L.  Croup  6D  -  Moving  Origins 

Up  to  3  moving  origins  may  be  defined.  Each  may  be  related  to  origins 
in  Groups  in,  6A  or  6B.  Moving  origins  may  be  read  in  from  tape  or 
generated  from  initial  position  and  rate  data.  Attitude  (pitch,  roll 
and  yaw)  data  of  ship  motion  may  be  read  in, 

M.  Group  7  -  Start  and  Stop  Times 

This  group  controls  the  time  spans  of  the  solution. 
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N.  Group  8  -  Refraction  Tables 

Up  to  10  refraction  index  tables  may  be  input.  A  ray  tracing  technique 
is  used  to  apply  tropospheric  and  ionospheric  refraction  corrections. 

O.  Group  9  -  Measurement  Definitions 

Up  to  40  measurements  (each  having  position,  position  arri  velocity, 
position,  velocity  and  acceleration)  of  the  types  mentioned  in 
Section  III  A  may  be  defined  and  base  line  time  delays  may  be 
defined. 

P.  Group  10  -  Measurement  Model 

Measurement  error  models  of  any  combination  of  those  listed  in 
Section  III  D  may  be  defined.  As  many  as  12  terms  may  be  defined 
for  each  of  the  40  measurements  provided  the  total  number  of 
estimable  terms  is  less  than  100. 

Q.  Group  11  -  Additional  Error 

Estimates  of  additional  random  error  may  be  defined. 

R.  Group  12  -  Measurement  Tape  Words 

The  locations  on  the  input  tape  or  input  cards  of  the  40  measurements 
and  their  random  errors  defined  in  Group  9  are  defined  in  this  group. 

S.  Group  13  -  System  Parameter  Estimates  for  this  Interval 

System  parameter  estimates  and  their  a  priori  estimates  of  error 
with  identification  numbers  less  than  101  may  be  estimated  or 
propagated.  Those  with  identification  numbers  larger  than  100 
are  not  estimated,  but  the  estimates  of  error  may  be  propagated. 

T.  Group  14  -  System  Control 

Control  of  measurement  entry  into  the  solution  is  available  via 
time  spans  and/or  minimum  elevation  cutoff.  Scaling  of  input 
measurements  is  also  available. 
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U.  Group  15  -  Tracking  Points 

Up  to  10  tracking  point  vector  pairs  may  be  defined. 

V.  Group  16  -  Initial  Approximations 

Up  to  5  sets  of  initial  approximations  may  be  defined  as  word  positions 
of  the  input  XYZ  or  AER  (position,  position  and  velocity  or  position 
and  velocity  and  acceleration)  or  right  ascension  and  declination. 

Or,  any  or  all  of  the  5  sets  may  be  single  points  input  on  cards. 

W.  Group  17  -  Tracking  Point  Orientations 

Up  to  20  sets  of  missile  attitude  information  may  be  input  to  be 
used  in  correcting  data  for  tracking  points. 

X.  Group  ISA  -  Free-Fall  Constants  the  Current  Interval 

Julian  day  number,  GMT  of  T-O,  free  fall  start  and  stop  times,  initial 
condition  data  and  integration  criteria  are  input. 

Y.  Group  18B  -  Solar  Flux  and  Geomagnetic  Index 
Tables  are  entered. 

Zl.  Group  18C  -  Drag  -  Lift  Information 

Tables  of  drag  and  lift  (up  to  300  entries  each)  may  be  entered,  along 
with  cross-sectional  area,  weight  and  reference  time  for  drag  and  lift 
error  model  polynomials.  The  tables  may  be  functions  of  time,  alti¬ 
tude  or  mach  number. 

Z2.  Group  18D  -  Weather  -  Atmosphere  Data 

The  built-in  atmosphere  profile  (density  and  sonic  speed  as  functions 
of  altitude  for  the  Patrick  standard  from  0  to  30  km,  and  U.S. 

Standard  1962  from  30  to  700  km)  may  be  called,  or  any  profile  of 
speed  of  sound,  density,  wind  speed  and  wind  velocity  may  be  entered. 
These  nay  be  entered  as  sets  polynomials  in  altitude  or  as  tables. 
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Z3.  Group  18E  -  Gravity  Model 

Built-in  geopotentials  may  be  called  or  any  geopotential  model  may 
be  read  in  on  cards.  Size  limits  are  order  24  degree  24  for  the  7004 
version  and  the  360/65  version,  and  order  30  degree  30  for  the  CYBER/ 
6500  version.  The  models  may  be  unnormalized,  APL  normalized  or 
Kaula  normalized.  The  models  may  be  truncated  at  specific  missile  or 
trajectory  altitudes. 

Z4.  Group  19  -  Plot  Scales 

Printer  plots  of  residuals  may  be  output  with  selectable  scales. 
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A.  Program  Certification 

Air  Force/RCA  policy  requires  that  each  production  computer  program 
(.any  program  that  processes  data  for  external  dissemination)  must  be 
certified,  fhe  certification  procedure  requires,  among  other  things, 
that  the  validity  of  the  program  results  be  verified  by  some  means 
external  to  the  program  such  as  nand  checkpoints  or  comparison  with 
a  certified  program.  All  production  versions  of  N1TE  are  certified. 

B.  Individual  Computer  Run  Validation 

rhe  validity  of  the  results  of  each  NITE  run  is  verified  by  evaluation 
of  the  following: 

1.  Convergence. 

Obviously,  the  results  are  of  no  use  unless  the  least  squares 
iteration  process  has  converged.  The  "outer  iteration"  is  the 
name  given  to  that  part  of  the  weighted  least  squares  process 
that  estimates  the  systematic  error  model  terms  (Groups  4B,  5A, 

13  and  the  initial  conditions  in  18A).  The  results  of  each 
iteration  are  printed,  i.e.,  the  corrections  applied  to  the 
estimates  of  the  previous  iteration,  and  the  estimates  of  the 
current  iteration.  Tne  analyst  in  charge  of  the  run  reviews 
the  magnitude  of  the  corrections  and  finds  the  run  acceptable 
if  the  magnitudes  of  the  corrections  become  insignificant  relative 
to  the  magnitude  of  the  current  estimate  and  its  uncertainty. 

2.  The  analyst  reviews  the  magnitude  of  the  final  systematic  error 
estimates  relative  to  the  a  priori  and  computed  systematic  error 
estimate  uncertainties.  Outliers  are  singled  out  for  further 
analysis  to  justify  their  values.  Usually  the  computed  uncertanties 
of  the  estimates  is  significantly  smaller  than  the  a  priori 
uncertainties,  a  reflection  of  the  strength  of  the  solution. 
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When  the  results  show  small  or  no  reductions  in  uncertainties, 
the  run  requires  more  analysis  to  verify  the  results  of  these 
"weak"  solutions. 

3.  The  trajectory  is  the  result  of  a  weighted  least  squares  combin¬ 
ation  of  the  data  corrected  for  systematic  error  model  terms  in 
what  is  called  the  "inner"  iteration  loop.  The  values  of  the 
computed  trajectory  are  analyzed  for  reasonableness  in  the  sense 
of  noise  content  and,  if  necessary,  comparisons  with  other  esti¬ 
mates  of  the  trajectory  (as  computed  from  other  single  sensors, 
or  other  NITE  runs). 

4.  The  correlation  matrix  of  the  systematic  error  model  terms  is 
reviewed.  Values  greater  than  1  are  indicators  of  matrix  inversion 
errors  which  can  happen  if  the  matrix  is  ill  conditioned  in  a  weak 
solution.  High  correlations  (less  than  or  equal  to  one)  are 
acceptable,  but  are  indicators  of  poor  separability  among 
coefficients. 

5.  The  magnitudes  of  the  residuals  are  reviewed.  The  residual  analysis 
outputs  means,  standard  deviations,  maximums  and  minimums  of  the 
residuals.  Means  of  the  residuals  significantly  different  than 
zero  require  analysis  for  the  cause: 

a)  A  priori  uncertainties  of  the  error  model  terms  may  be  small 
enough  to  give  the  a  priori  value  too  much  weight,  preventing 
full  resolution  of  the  systematic  error.  Action:  raise 

a  priori  uncertainty  or  change  c’ne  estimate  of  the  systematic 
error  term. 

b)  The  dynamic  weighting  may  be  weighting  out  trended  residuals. 
Action:  verify  weighting. 

c)  Ihe  systematic  error  model  may  be  incorrect.  Action:  review 
trends  in  residuals,  find  justification  for  new  model  from 
analysis  of  system  performance. 


23 


Review  the  Distribution  of  the  Residuals 

This  is  best  accomplished  by  review  of  the  plots  of  the  residuals. 
Ideally,  the  residual  distribution  should  be  random  about  zero. 

Any  deviations  should  be  justified  by  analysis  of  the  weighting 
(high  random  error  estimates  hence  low  weights  indicate  marginal 
system  performance,  low  elevation  tracking  being  contaminated  by 
refraction  errors,  multi-path,  clutter,  etc.),  and  a  review  of 
the  modeling  . 


VI.  A  PRIORI  INFORMATION 


A.  Estimates  of  the  uncertainty  of  error  model  coefficients  are  derived 
primarily  by  the  Technical  Analysis  group.  Historical  statistical 
histories  are  maintained  of  the  results  of  all  BET  and  satellite 
solutions.  The  error  model  coefficients  and  their  uncertainties 

are  evaluated,  surmarized  and  the  results  are  produced  in  a  "Quarterly 
Accuracy  Report"  (QAR).  Uncertainties  of  any  terms  not  presented  by  tie 
report  are  derived  by  the  analyst  in  charge  and  are  based  upon  the 
facts  available,  such  as  the  expected  magnitude  of  the  coefficient, 
the  strength  of  the  solution,  the  quality  of  the  measurement  and  the 
information  presented  by  NITE  in  the  partial  derivative  output. 

B.  Random  Error  Estimates 

Estimates  of  random  error  used  in  weighting  are  usually  derived  from 
scaled  differences  between  smoothed  and  unsmoothed  measurements.  The 
smoothing  span  of  the  filter  is  selected  to  ensure  maximum  smoothing 
without  the  removal  of  missile  motion.  The  scaled  differences  are 
augmented  by  as  much  as  57a  of  the  refraction  correction,  to  ensure 
that  range  and  elevation  are  properly  weighted  to  account  for  residual 
refraction,  multipath  and  ground  clutter  errors.  The  random  error 
estimates  are  monitored  to  ensure  that  they  are  within  the  range  of 
past  performance  of  that  system  as  indicated  in  the  Tech  Analysis  QAR. 
Occasionally,  a  system  is  used  that  cannot  be  smoothed.  Q^R  estimates 
are  used  if  available;  if  not,  estimates  are  derived  from  comparisons 
with  more  accurate  systems,  or  by  analysis  of  NITE  residuals  where 
the  system  in  question  is  weighted  out  (provided,  of  course,  that  the 
BET  is  known  to  have  less  noise  than  the  system  in  question).  If 
data  are  collected  while  the  tracked  object  is  in  free  fall,  and  free- 
fall  constraints  are  used  in  NITE,  a  good  estimate  of  random  noise 
can  be  found  in  the  residuals. 
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Robert  W.  Mai 

US  Army  Yuma  Proving  Ground,  Yuma,  Arizona 


ABSTRACT 

This  paper  describes  the  background  for  the  development  of  US  Army 
Yuma  Proving  Ground  (USAYPG)  Best  Estimate  of  Trajectory  (BET)  program. 

The  program  is  called  GPSBET  which  employs  optimal  estimation  techniques 
to  combine  measurements  from  several  sources  and  arrive  at  a  single 
trajectory  estimate.  The  estimation  technique  employed  is  Kalman 
filtering  and  smoothing.  The  implementation  of  the  filter  and 
smoothing  equations  is  described  in  this  paper,  as  well  as  the  ancillary 
processing  performed  to  prepare  the  measurements  for  filtering,  '"he 
input  to  and  output  from  GPSBET,  in  addition  to  the  validation  proced'^es , 
are  also  described.  The  performance  of  GPSBET  in  terms  of  accuracy  and 
computer  run  time  is  also  documented, 
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1.  HISTORY 

In  the  pre-laser  era,  YP6  used  a  multistation  cinetheodol ite  system 
to  obtain  data  for  a  Best  Estimate  of  Trajectory.  The  trajectory  was 
obtained  from  an  odle  solution  (minimizing  the  sum  of  squares  of  dis¬ 
tance  residuals),  followed  by  a  Davis  solution  (minimizing  the  sum  of 
squares  of  angular  residuals),  followed  by  a  second  degree  moving  arc 
smoothing  routine.  In  1973,  YPG  received  a  PATS  (Precision  Aircraft 
Tracking  System)  Laser.  Upon  comparison  of  laser  derived  trajectories 
with  cinetheodol ite  derived  trajectories,  it  became  apparent  that  the 
laser  had  the  potential  for  providing  position  data  In  a  very  automated 
mode  and  with  accuracy  as  good  as  cinetheodol ites  but  dependent  only  on 
distance  from  the  PATS  rather  than  being  complicated  by  geometry.  A 
software  development  effort  was  initiated  in  the  fall  of  1973  to  comoute 
trajectory  using  some  combination  of  laser  and  cinetheodol ite  data. 

The  software  development  was  based  on  Kalman  filtering  of  measurements. 
The  Kalman  filter  was  ideal  for  ease  of  accepting  any  combination  of 
range  or  azimuth  or  elevation  measurements  subject  only  to  the  restriction 
that  they  come  in  chronological  order.  Subsequent  analysis  of  trajectory 
data  and  measurement  residuals  revealed  that  the  range  measurements 
from  the  laser  were  the  strength  of  any  laser/cine  solution.  YPG  then 
purchased  two  additional  lasers.  Further  tests  and  simulations 
demonstrated  that  cinetheodol ite  data  contributes  little  or  nothing  to 
the  accuracy  of  a  3-laser  solution  and  was  very  expensive  to  acquire 
and  reduce  and  it  slowed  data  turnaround  tremendously. 

In  mid-1974,  the  Aerospace  Corporation,  in  support  of  the  NAVSTAR 
Global  Positioning  System,  helped  YPG  with  the  develooment  of  Kalman 
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filter  equations  including  the  incorporation  of  INS  (Inertial  Navigation 
System)  and  IMU  (Inertial  Measurement  Unit)  data.  INS  data  consists  of 
sensed  velocities  in  three  components  and  IMU  data  consists  of  gimbal 
angles  which  reflect  the  pitch,  roll  and  yaw  of  an  aircraft.  The  purpose 
of  INS  measurements  was  to  lower  the  uncertainty  in  velocity  estimates. 
The  purpose  of  the  IMU  measurements  was  to  improve  lever  arm  corrections 
by  providing  accurate  attitude  Information . 

2.  TERMINOLOGY 


YPG's  Best  Estimate  of  Trajectory  program  is  called  GPSBET.  The 
estimation  is  accomplished  using  Kalrr.an  filter  equations  and  the  Rauch- 
Tung-Stnebel  recursive  smoothing  equations.  It  is  appropriate  therefore 
to  define  the  components  of  the  state  vector  which  are  estimated  by  the 
Kalman  filter/smoother.  The  state  vector  consists  of  X,  Y,  Z  components 
of  position,  velocity  and  acceleration  plus  a  bias  state  for  each 


measurement  type  from  each  instrument.  For  a  3-laser  solution,  the 


state  vector  has  18  components: 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 


X 

Y 
Z 
X 

Y 
Z 
X 

Y 
Z 


Range  bias 
Azimuth  bias 
Elevation  bias 
Range  bias 
Azimuth  bias 
Elevation  bias 
Range  bias 
Azimuth  bias 
Elevation  bias 


Position 
Position 
Position 
Velocity 
Velocity 
Velocity 
Acceleration 
Acceleration 
Acceleration 
Laser  1 
Laser  1 
Laser  1 
Laser  2 
Laser  2 
Laser  2 
Laser  3 
Laser  3 
Laser  3 


When  INS/ IMU  measurements  are  used  then  6  more  bias  states  are  added: 
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19 

Orientation  bias 

East  axis 

20 

Orientation  bias 

North  axis 

21 

Orientation  bias 

Vertical  axis 

22 

Velocity  bias 

East  axis 

23 

Velocity  bias 

North  axis 

23 

Velocity  bias 

Vertical  axis 

The  state  vector  estimates  of  position,  velocity  and  acceleration 
are  in  a  right  handed  cartesian  coordinate  system  with  origin  at  a 
specific  point  on  the  ground  at  YPG.  The  bias  estimates  for  each  laser 
are  in  a  spherical  coordinate  system  centered  at  the  respective  lasers 
with  zero  azimuth  being  True  North  at  each  laser  and  zero  elevation 
being  horizontal  at  each  laser.  The  orientation  biases  are  relative  to 
the  gyro  wander  angle  reference  frame  as  are  the  velocity  bias  estimates. 
3.  ANCILLARY  PROCESSING 

The  purpose  of  trajectory  estimation  at  YPG  is  generally  to  provide 
a  position,  velocity  and  acceleration  time  history  for  a  particular  point 
on  a  particular  target  in  a  standard  coordinate  system.  Therefore,  a 
trajectory  estimation  program  must  include  modules  for  atmospheric  cor¬ 
rections,  systematic  error  correction,  coordinate  conversion,  and  lever 
arm  corrections,  in  addition  to  the  noise  reduction  and  state  vector 
estimation  algorithms.  Also,  the  actual  measurements  going  into  a  BET 
program  are  never  perfect  and  the  imperfections  are  not  time  stationary 
so  that  a  considerable  amount  of  checking  for  wild  points  and  for  changes 
in  noise  distribution  must  be  done.  Sections  3.1  through  3.7  discuss 
some  of  the  special  data  handling  and  evaluation  procedures  that  are 
included  in  GPSBET  to  help  ensure  that  the  Kalman  filter  parameters 
adequately  describe  the  data  being  filtered. 

3.1  ATMOSPHERIC  CORRECTIONS 

GPSBET  corrects  all  laser  range  and  elevation  measurements  for 
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atmospheric  refraction.  The  model  currently  being  used  is  a  two  param¬ 
eter  model  where  the  parameters  are  computed  from  RAWINDSONDE  data 
obtained  generally  within  two  hours  of  the  time  for  which  a  trajectory 
estimate' is  desired.  The  model  assumes  a  constant  temperature  lapse 
rate  so  that  when  a  target  is  located  in  or  near  a  low  level  temperature 
inversion,  significant  errors  can  result  (where  significant  means  any¬ 
thing  over  .5  meter  ) . 

3.2  TRACKING  SYSTEM  ERRORS 

The  laser  measurements  are  corrected  for: 

Range  bias  - 

Range  scale  factor  -  Rs 

Azimuth  bias  -  AP 

Elevation  bias  -  Ek 

North  tilt  -  NT 

East  tilt  -  £T 

Collimation  -  CJ5L 

These  corrections  are  made  by  computing  coefficients  during  calibration. 
Calibration  includes  collecting  at  least  400  measurements  from  each 
laser  as  each  lases  on  eight  targets  at  45°  increments  around  each  laser. 

The  equations  for  correcting  measurements  (Rm,  Am,  Em)  to 
(Rc.  Ac,  Ec)  are: 

03  ^c  "  Rm  +  ^s  *  Rm  +  Rb 

Ec  =  Em  +  Eb 

03  Ac  =  Am  +  Ab  +  sin'1  |jin  jCOL) j 

C3l  Ec  =  sin'1  {sin  (Ec)  *  cos  (COL)j 

Mislevel  of  each  laser  is  accounted  for  by  rotation  of  the  predicted 
state  X,  Y,  Z  from  the  common  cartesian  coordinate  system  to  the  laser 
misleveled  system  (which  includes  earth's  curvature  from  the  selected 
origin  to  the  laser  sites).  This  is  a  convenient  and  precise  method  of 
obtaining  azimuth  and  elevation  predictions,  Ap  and  Ep,  which  are  in  the 


36 


same  reference  system  as  the  Ac  and  Ec  computed  in  QQ  and  £33. 

Other  sources  of  tracking  system  errors  have  been  investigated. 

Encoder  eccentricity  errors  are  observable  and  although  they  are  small, 
will  be  included  In  the  BET  when  the  calibration  software  has  been 
modified  to  compute  coefficients.  Lens  sag  is  negll gable  for  the 
YPG  lasers.  Servo  lag  is  negligable  as  long  as  mount  acceleration  is 

p 

significantly  less  than  80  mill iradians/sec  (which  is  generally  the 
case  for  aircraft  tracking).  GPSBET  does  estimate  range,  azimuth  and 
elevation  rates  and  accelerations  which  can  be  used  to  model  servo 
response.  The  model  appropriate  for  YPG's  lasers  is: 

L  =  K  *  A, 

where 

L  =  Angular  lag 

K  =  Constant  dependent  upon  qa in  of  servo  system 

A  =  Angular  acceleration  experienced  by  the  mount 
3.3  EARTH'S  CURVATURE  CORRECTIONS 

GPSBET  produces  state  vector  estimates  in  a  cartesian  coordinate 
system  with  origin  at  input  latitude,  longitude  and  altitude.  The 
geodetic  coordinates  of  the  origin,  as  well  as  the  tracking  instruments, 
are  input  relative  to  a  reference  ellipsoid  whose  parameters  can  be  varied. 
The  model  of  the  earth  used  is  Clarke's  Spheroid  of  1866  with  the  WGS 
72  Datum.  At  each  measurement  time,  the  estimate  for  the  state  vector 
at  that  time  is  corrected  for  earth's  curvature  by  rotating  the  X,  Y,  Z 
through  the  difference  in  latitude  and  longitude  from  the  input  origin 
to  tne  instrument  location  (plus  the  mislevel  of  the  instrument).  The 
rotated  X,  Y,  Z  can  then  be  transformed  to  spherical  coordinates  (range, 
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azimuth,  elevation)  and  compared  with  the  measurements  (which  have  been 
corrected  for  tracking  system  errors  as  described  in  3.2).  The  Kalman 
filter  uses  the  difference  between  predictions  and  measurements  (called 
residuals)  to  update  state  vector  and  covariance  matrix  estimates  and 
to  make  edit  checks. 

3.4  ATTITUDE  ESTIMATION  AND  LEVER  ARM  CORRECTIONS 

Range  instrumentation  can  not  usually  provide  measurements  of 
the  position  or  velocity  of  the  particular  point  on  a  target  for  which 
a  reference  trajectory  is  required.  In  the  case  of  laser  tracxers,  a 
retrorefl ector  package  must  be  installed  on  the  target  and  the  location 
of  that  package  cannot  coincide  with  the  desired  reference  point.  For 
radar,  the  same  situation  exists  for  beacon  tracking.  In  the  case  of 
cinetheodol ite  tracking ,  the  desired  reference  point  is  not  always 
identifiable  on  film.  When  very  accurate  position  and  velocity  informa¬ 
tion  is  required  for  a  particluar  point  on  a  target  (such  as  an  antenna 
or  the  center  of  gimbals  in  a  gyro  package),  then  corrections  must  be 
made  to  allow  for  the  separation  of  the  instrumentation  reference  point 
and  the  desired  reference  point.  On  aircraft,  the  separation  may  be 
on  the  order  of  10  meters.  When  .5  meter  accuracy  is  required,  then 
10  meter  separation  of  reference  points  must  be  carefully  removed.  This 
is  done  in  GPSBET  by  continual ly  updating  estimates  of  the  attitude 
matrix.  The  attitude  matrix  is  a  rotation  matrix  to  go  from  the  aircraft 
body  coordinate  system  to  the  reference  cartesian  coordinate  system. 

It  includes  pitch,  roll,  yaw  information  about  the  aircraft  plus  earth's 
curvature  from  the  aircraft  to  the  origin  of  the  reference  cartesian 
system. 


To  accomplish  lever  arm  corrections,  the  following  vectors  and 


matrices  must  be  computed  or  input: 


M  =  Rotation  matrix  to  correct  for  Earth's  curvature  from  the 

origin  of  the  reference  cartesian  system  to  a  cartesian  system 
with  aircraft  as  origin. 

=  Matrix  transpose  of  M. 


PRY  =  Rotation  matrix  to  correct  from  pitched,  rolled  and  yawed 
aircraft  body  coordinate  system  to  locally  horizontal 
XY-plane  (Z  up)  at  the  aircraft. 


Attitude  matrix  to  express  rotation  from  aircraft  body 
coordinate  system  to  reference  cartesian  system. 

Vector  containing  distances  in  aircraft  body  coordinate  system 
from  desired  reference  point  to  instrumentation  reference 
point. 


Kalman  filter  estimates  of  position  of  desired  reference  point 
on  aircraft. 


Position  of  instrumentation  reference  point  on  the  aircraft. 


To  incorporate  a  measurement  from  a  particular  instrument,  6PSBET  does 


the  following  for  lever  arm  correction: 


where 


and 


M 


M 

r x 

n 

=  A  *  I  Y 

[ ZI 

1  Z 

L  J 

A 

=  MT  *  PRY 

cDL  sDLsPH 

sPHIsDL  cPHlcPH  +  sPHlcDLsPH 

cPHlsDL  sPHIcPH  -  cPH IcDLsPH 


-sDLcPH 

c PH  I s PH  -  sPHIcDLcPH 
sPHIsPH  +  cPHlcDLcPH 
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cBcP 

-cBsPsR  +  sBcR 

-cBsPcR 

-  sBsR 

PRY  = 

-sBcP 

sBsPsR  +  cBcR 

sBsPcR  - 

cBsR 

sP 

cPsR 

cPcR 

sXXX  -  sin  of  the  angle  XXX 
cXXX  -  cos  of  the  angle  XXX 

DL  -  difference  in  longitude  of  aircraft  and  origin 
PH  -  Latitude  of  origin 
PHI  -  latitude  of  aircraft 
P  -  pitch  of  aircraft 
R  -  roll  of  aircraft 

B  -  heading  of  aircraft  with  respect  to  True  North 
The  pitch,  roll  and  heading  of  the  aircraft  are  obtained  from  gimbal 
angle  measurements  when  available.  Otherwise,  velocity  and  acceleration 
estimates  are  used  to  estimate  pitch,  roll  and  heading. 

The  attitude  matrix.  A,  computed  at  each  measurement  time  is  an 
important  part  of  the  output  of  GPSBET.  It  can  be  used  to  determine  the 
position,  velocity,  acceleration  of  a  point  on  the  aircraft  other  than 
that  estimated  by  the  state  vector. 

3.5  RETROREFLECTOR  IDENTIFICATION 

The  NAVSTAR  Global  Positioning  System  project  is  currently  the 
principal  user  of  YPG's  laser  tracking  instrumentation  and  the  GPSBET 
program.  For  GPS  testing,  two  retroref lector  packages  are  mounted  on  some 
aircraft  (e.g.  ,  C-141  and  F-4).  The  purpose  is  to  eliminate  blockage  of 
retrorefl ectors  by  the  wings,  fuselage,  etc.  This  results  in  virtually 
complete  coverage  from  three  laser  trackers  throughout  all  maneuvers. 

Since  the  retros  are  separated  by  5-10  meters  from  each  other,  the 


GPSBET  program  needs  an  algorithm  to  decide  which  retro  is  being  tracked 
by  a  given  laser  so  that  the  appropriate  lever  arm  corrections  might  be 
applied.  The  alogrithm  being  employed  involves  computing  the  look  angle 
from  a  given  laser  to  the  axis  of  each  retroreflector  package.  The  retro 
packages  used  are  hemispherical  and  have  been  mounted  on  the  top,  pointing 
up,  and  the  bottom,  pointing  down,  of  an  aircraft.  The  retro  whose  axis 
makes  the  smallest  angle  with  the  line  of  sight  to  a  given  laser  is  the 
retro  Identified  as  the  most  likely  source  for  returning  the  light  to  the 
laser  and  lever  arm  corrections  are  made  from  the  reference  point  chosen 
for  the  state  vector  to  the  selected  retro. 

The  use  of  two  retroreflector  packages  does  solve  the  blockage 
problem  (i.e.,  all  lasers  almost  always  have  line  of  sight  to  one  retro 
or  the  other^but  causes  degradation  of  accuracy  since  the  identification 
is  imperfect  and  lasers  can  get  light  returned  from  both  retro  packages 
at  the  same  time.  The  identification  scheme  is  very  sensitive  to  errors 
in  the  estimated  roll  of  the  aircraft  and  this  gives  rise  to  the  erroneous 
identification  of  retro  switches  which  in  turn  puts  noise  on  velocity 
and  acceleration  estimates. 

3.6  EDITING  AND  REINITIALIZATION 

The  GPSBET  program  accomplishes  editing  by  checking  the  magnitude 
of  measurement  residuals  prior  to  the  incorporation  of  the  measurement 
into  the  filter.  The  procedure  includes  the  following  steps: 

1.  Compute  the  difference  between  the  measurement  and  the  filter's 
prediction  for  the  measurement  (this  difference  is  called  the  measure¬ 
ment  residual  and  is  computed  as  described  in  3.3). 

2.  Compute  the  expected  variance  of  the  residual  by  adding  the 
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variance  in  the  prediction  (computed  from  the  covariance  matrix  of  the 
state  vector)  to  the  variance  of  the  measurement  noise  (a  priori 
information  required  as  input  to  GPSBET). 

3.  Divide  the  residual  by  the  uncertainty  in  the  residual  (the 
square  root  of  the  variance  described  in  2.). 

4.  If  the  residual  is  more  than  N  times  the  uncertainty  in  the 
residual  (where  N  is  an  input  constam  _ually  set  to  5)  then  the 
measurement  is  not  included,  otherwise  the  Kalman  filter  equations  are 
executed  to  correct  the  predicted  state  vector  and  covariance  matrix 
estimates. 

5.  If  editing  has  occurred,  then  counters  are  incremented  for: 

-  total  numoer  of  edits  of  measurements  of  this  type 

-  number  of  consecutive  edits  of  measurements  of  this  type 

-  number  of  consecutive  edits  of  measurements  from  this  instrument. 

The  edit  counters  described  in  5.  are  checked  throughout  the 

filtering  of  data  and  when  input  l  imns  are  exceeded,  some  action  is 
taken,  e.g., 

-  the  filter  is  reinitialized  when  at  least  two  of  the  three  types 
of  measurements  from  all  lasers  have  oeen  edited  for  an  input  number  of 
seconds 

-  the  filter  is  reinitialized  when  at  least  two  of  the  three 
measurement  types  have  been  edited  for  an  input  number  of  consecutive 
measurement  times  for  a  given  laser 

-  the  filter  is  reinitialized  when  a  given  measurement  is  edited 
consecutively  more  often  than  an  input  limit 

-  printing  of  edit  messages  is  stopped  when  an  input  limit  is  reached. 


Reinitialization  is  accomolished  by  converting  the  range,  azimuth, 

elevation  data  to  X,  Y,  Z  in  the  common  cartesian  system  for  each  of 

the  lasers.  When  three  consecutive  second  differences  are  within  an 

input  criterion  for  X,  Y  and  Z  of  a  given  laser,  then  a  second  degree 

fit  is  made  to  those  five  X  coordinates,  five  Y  coordinates  and  five  Z 

coordinates.  The  polynomial  and  its  derivatives  are  evaluated  at  the 

time  of  the  first  measurement  in  the  fit.  This  provides  X,  Y,  Z,  X,  Y, 

•  •  •  •  •  •  • 

Z,  X,  Y,  and  Z  which  are  used  as  initial  state  vector  components  for 
the  Kalman  filter.  The  covariance  matrix  is  reinitialized  to  the  same 
values  that  were  input  at  the  beginning  of  the  estimation.  This  scheme 
for  reinitial izing  the  filter  is  usually  used  to  get  initial  state 
vector  estimates  at  the  beginning  of  the  trajectory  estimation. 

An  eleven-point  fit  is  an  option  for  initialization  and  reini¬ 
tialization,  but  is  rarely  required;  since  the  Kalman  filter  initial¬ 
ization  transients  dampen  before  highlv  accurate  estimates  are  required. 
In  most  cases;  the  initial  uncertainties  in  the  state  vector  are  input 
as  large  numbers  (e.g., 1000  meters,  500  meter/sec,  100  meters/sec/ 
sec  for  position,  velocity,  acceleration  components)  so  that  the 
Kalman  filter  will  accept  measurements  even  though  the  initialization 
was  not  very  good. 

3.7  RESIDUALS  SUMMARIES 

Probably  the  most  important  output  from  the  GPSBET  program  from 
an  analytical  standpoint,  (i.e.,the  output  that  provides  the  best 
measure  of  filter  performance,  measurement  quality  and  quantity  and 
accuracy  of  a  priori  information)  are  residuals  sunmary  reports. 


These  summaries  are  output  at  a  rate  determined  by  input  parameters 
and  can  be  selected  to  come  at  a  particular  time  interval  or  after 
a  certain  number  of  measurements  have  been  received.  After  a  set  of 
measurements  has  been  processed  for  a  given  instrument  (e.g.,  range, 
azimuth,  elevation  from  laser  1 ),  the  estimates  for  X,Y,Z  are  translated 
to  the  appropriate  instrument,  corrected  for  earth's  curvature, 
corrected  to  the  instruments  misleveleo  state,  converted  to  spherical 
coordinates  and  differenced  with  the  calibrated  measurements  to  form 
residuals.  Only  for  those  measurements  not  edited  are  the  residuals 
tabulated.  The  tabulated  residuals  are  summed,  their  squares  are 
summed  and  the  number  of  unedited  measurements  of  each  type  is  counted. 
At  requested  intervals  the  mean  of  the  residuals,  RMS  of  the  residuals 
and  number  of  measurements  used  since  the  last  report  is  printed 
for  each  measurement  type  for  each  instrument.  The  number  of  measure¬ 
ments  used  since  the  last  report  tells  how  much  editing  and  data 
dropouts  were  occurring.  The  RMS  of  residuals  is  an  excellent  measure 
of  the  noise  on  the  measurements  and  is  used  to  validate  the  a  priori 
information.  The  mean  of  residuals  should  be  near  zero,  but  when  it 
isn’t  points  to  problems  such  as  erroneous  calibrations,  misidenti- 
fication  of  retros,  lack  of  sufficient  filter  response  to  process 
noise,  inconsistent  measurement  sources,  etc. 

4.  KALMAN  FILTER 

The  program  GPSBET  uses  Kalman  filter  equations  to  derive  a 
trajectory  estimate  at  each  time  for  which  a  measurement  is  made 
available.  Section  4.1  explains  the  notational  conventions  to  be 
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used  in  4.2  to  describe  the  input  parameters  controlling  the  behavior 
of  the  Kalman  filter  and  in  4.3  to  outline  the  implementation  of  the 
Kalman  filter  equations  in  GPSBET. 

4.1  NOTATION 

Capitol  letters  (e.g.,A,M,X)  will  be  used  to  represent  matrices. 
Sometimes  these  matrices  have  only  one  column  or  one  row,  in  which  case 
they  will  be  referred  to  as  vectors, but  no  notational  distinction  will  be 
made.  Time  will  be  denoted  by  t,  and  tn  refers  to  the  nth  time  in  a 
sequence.  At  will  refer  to  a  time  interval  from  some  time  tm  to 
another  time  tn  where  the  order  of  the  subtraction  will  be  made  clear 
in  the  text.  When  a  capital  letter  is  subscripted  such  as  Xn  this 
will  refer  to  the  estimate  of  X  at  time  tn.  A  vector  (matrix)  symbol 

A  a 

with  a  carrot  over  it  (e.g.,  Xn,Pn)  refers  to  a  predicted  value  for  the 
vector  (matrix)  at  time  tn.  In  all  equations, the  symbol  *  refers  to 
multiplication.  Lower-case  letters  (e.g.;  a,b,n,j)  refer  to  scalars. 

The  lower  case  letters  i,j,k,l,m,n  are  reserved  for  integer  scalars. 

The  scalars  x,y,z,x,y,£,x,y,z  are  reserved  for  the  position, 
velocity,  acceleration  components  in  a  particular  state  vector  estimate. 

The  notation  b*X  refers  to  multiplication  of  a  scalar  and  a 
matrix  (in  particular,  At  *  0  is  the  mul tinl ication  of  the  matrix 
Q  by  the  scalar  At  which  is  a  time  interval).  Multiplication, 
addition  and  subtraction  in  matrix  equations  implies  that  the  number 
of  rows  and  columns  in  each  matrix  allows  such  operations.  The 
notation  refers  to  the  transpose  of  the  matrix  H. 


Xn  -  Kalman  filter  estimate  for  the  state  vector  after  incorporation 
of  measurements  from  a  given  instrument  at  time  tn  (the  first  9 
components  are  x,y,z,x,y  ,z  ,x,y,z) 

Pn  -  Kalman  filter  estimate  for  the  covariance  matrix  of  the  state 
vector  Xn  (the  diagonal  of  Pn  contains  tne  estimates  of  the  variance 
of  the  components  of  the  state  vector) 

$n,n+l  '  Transformation  matrix  to  estimate  a  value  for  Xn+i  at  time 
tn+1  given  Xn  and  At=tn+i-tn 

A 

Xn^  -  Prediction  of  Xn+i  given  Xn  and  A  t*tn+i -tn 

/i 

*n+l  =  $n’n+l  *  ^n 

Q  -  Effect  of  process  noise  on  covariance  matrix  P  per  unit  of  time 

A 

Pn+1  -  Prediction  of  Pn+]  given  Pn  andZ)t  =  tn+-]  -  tn 
^n+1  =  Sn t  n+1  *  (^n  +^>t  *  Q)  Snt  n+l 

9n  -  Measurement  taken  at  time  tn 

A  A 

Qn  -  Prediction  of  0n  given  Xn 

An  -  Transformation  matrix  to  go  from  state  vector  domain  to  measure¬ 
ment  domain 

A  A 

®n  =  ^n  * 

NOTE:  In  GPSBET,  this  transformation  is  accomplished  with  the  trigono¬ 
metric/algebraic  relationships  between  cartesian  coordinates 
(X,  Y,  Z)  and  spherical  coordinates  (range,  azimuth,  elevation). 

The  transformation  of  state  vector  uncertainties  to  uncertainties 

A  T 

in  9n  are  estimated  by  Hn  *  Pn  *  Hp  where  Hn  is  described  below. 

Hn  -  Partial  derivatives  of  measurement  to  be  incorporated  at  time 
tn  (i.e.,  9n)  with  respect  to  the  state  vector: 
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Kn  -  Gain  vector  for  Kalman  filter  computed  for  incorporating  the 
measurement  0n  taken  at  time  tn 

r  -  Scalar  denoting  the  measurement  residual  i.e., 

r  3  ®n  -  9n 

4.2  KALMAN  FILTER  INPUT  PARAMETERS 

This  section  will  describe  only  those  input  parameters  that  are 
required  in  the  Kalman  filter  implementation  in  GPSBET.  Many  other  param¬ 
eters  are  input  to  the  program  for  controlling  data  flow,  manipulating 
input  data,  selecting  processing  options,  selecting  output  options,  etc. 

These  other  parameters  will  be  described  in  Appendix  A  so  as  not  to 
clutter  the  main  body  of  this  report. 

D  -  An  array  containing  the  standard  deviations  of  the  noise  on 
each  measurement  type  for  each  instrument.  It  is  double  dimensioned  so 
that  0  (i,  j)  is  the  standard  deviation  of  the  i^  measurement  type  of  the 
jtU  instrument.  For  example,  D  (2,  3)  -  .000200,  could  mean  the  standard 
deviation  of  the  azimuth  measurements  for  the  third  laser  is  200  microradians . 
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Q  -  A  matrix  the  same  size  as  the  covariance  matrix  P, but  which 
contains  the  increased  variance  expected  in  the  state  vector  due  to 
process  noise  expected  in  1  second. 


P0  -  Initial  values  of  tne  covariance  matrix,  usually  set  very  large 
so  that  data  is  accepted  and  only  consistent  data  will  drive  the  diagonal 
of  P  to  small  values. 

XQ  -  Initial  values  tor  the  state  vector  which  can  be  input  or 
computed  from  laser  data  as  oescrioed  in  Section  3.6. 

4.3  FILTER  EQUATIONS 

To  describe  the  implementation  of  the  Kalman  filter  equations  in 
GPS8ET,  this  section  will  start  with  a  state  vector  estimate  Xn  anG  an 
estimate  of  its  covariance  matrix  Pn  at  time  tn.  The  equations  and 
intermediate  steps  performed  by  GPSBtT  to  incorporate  a  measurement 
9n+1  taken  at  time  tn+i  will  be  described.  The  estimates  Xn  and  Pn  may 
have  been  initial  estimates  (denoted  XQ  and  PQ)  as  described  in  Section 
4.2.  The  measurement  0n+j  is  assumed  to  have  been  corrected  for  atmos¬ 
phere  effects  as  described  in  Section  3.1  and  for  tracking  system  errors 
as  described  in  Section  4.2. 

Step  .  The  estimates  for  the  state  vector,  Xn,  and  covariance  matrix, 
Pn,  at  time  t^,  are  projected  ahead  (or  back  if  measurements  come  in 
reverse  chronological  order)  to  obtain  predictions  Xn+i,  Pn+j  for  their 
value  at  time  tn+1 : 

A*  *  Vl  *  ln 

An+)  =  Sn*  n+1  *  Xn 

^♦1  *  V  n+1  (pn  +  At  *  Q)  *  sj.  n+1 
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t2 

00  t  0  0  0  0  0  .  . 

10  0  t  0  0  0  0  .  . 

0  10  0  t  0  0  |  0  .  . 

00100  t  0  0  0  .  . 

000100  too.. 
0000100  to.. 
000001000.  .  . 


Step  2:  The  predicted  value  for  the  state  vector,  Xn+i,  is  transformed 

A 

to  a  predicted  value,  9n+i,  for  the  measurement,  9n+l-  The  computation 

A 

of  9n+i  is  accomplished  with  algebra/trigonometry  rather  than  matrix 
multiplication.  For  example,  if  9n+j  is  an  azimuth  measurement  from 

A 

laser  1,  then  the  X,  Y,  Z  components  of  Xn+i  are  rotated  and  translated 

to  the  local  cartesian  system  at  laser  1  (call  them  X'  ,  Y'  ,  Z'  ) .  Then 
A  A  .  • 

9n+i  is  computed:  9n+i  =  sin  1  (  ,■  -~==:  A  )  +  b, where  b  is  the 

A  ^  (x‘  )  +  (y‘  )2 

component  of  Xn+j  that  represents  the  filter's  estimate  of  azimuth  bias 
for  laser  1.  The  measurement  bias  estimates  remain  zero  unless  the 
initial  uncertainty  in  the  measurement  bias  states  is  large  enough  and 
measurements  from  other  sources  are  good  enough  to  allow  the  Kalman 
filter  gain  for  measurement  bias  to  cause  non-zero  biases  to  be  estimated. 
Since  YPG's  laser  trackers  are  calibrated  for  every  mission  with  sur¬ 
veyed  targets, the  GPSBFT  program  is  usually  forced  not  to  estimate  laser 
measurement  biases  by  setting  the  initial  uncertainty  very  small.  For 
INS  velocity  measurements  on  the  other  hand  biases  of  1.-2.  meters/ 
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second  are  usually  estimated.  These  biases  result  from  gvro  drift  and 
are  modeled  rather  carefully  using  the  gyro  specifications  and  knowledge 
of  the  dynamics  experienced  by  the  gyro. 

T  A  A 

The  above  procedure  for  computing  9n+1  from  Xn+1  will  be  denoted 
subsequently  by 

®n+l  =  ^n+1  *  Xn+j  , 

although, as  described  above, a  matrix  An+i  does  not  have  to  be  found;  since 
algebraic  and  trigonometric  relationships  can  be  used  instead,  and  tnese 
relationships  provide  an  exact  transformation. 

Step  3:  The  difference  between  the  measurement  0n+i  and  the  prediction 

A 

9n+2  is  computed: 

A 

r  =  ®n+l  "  ®n+l 

The  scalar  r  is  called  the  residual. 

Step  4 :  The  standard  deviation  of  r,  denote  it  by  sr,  is  computed  next. 
Since  r  is  the  difference  of  two  quantities  presumed  to  be  statistically 
independent  (reasonable  since  the  measurement  at  time  tn+i  has  not  been 
used  at  all  to  arrive  at  9n+i),  then 

4  »  d2  +  [d  (i,  j)]2 

where 

D  (i,  j)  -  a  priori  estimate  of  the  standard  deviation  of  measure¬ 
ment  type  i  from  instrument  j,  and  9n+i  is  a  type  i  measurement  from 
instrument  j. 

A 

d  -  uncertainty  in  the  estimate  9n+i,  which  can  be  computed  from 
the  uncertainty  in  Xn+j  and  the  partial  derivatives  of  measurement  9n+j 

A 

with  respect  to  Xn+j : 
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d2  1  ”n+l  *  ?n+l  *  hJ+1 


Thus 


Step  5:  An  edit  check  is  made  next.  If 


< 


l  , 


£  is  input) 


then  the  remaining  Kalman  filter  eouations  are  executed.  If 


?  > 

V  — 


4, 


£  is  input) 


then  the  measurement  9n+i  is  edited.  The  limit  on  the  number  of  standard 
deviations  allowed  (denoted  by  £  above)  is  an  input  parameter  (often 
set  to  5.).  When  editing  is  indicated  all  appropriate  edit  counters  are 
incremented,  an  edit  message  may  be  printed,  the  following  assignments 
are  made: 


A 

*n+I  =  *n+l 
^n+1  =  ^n+1 

and  the  Kalman  filtering  cycle  for  measurement  9n+i  is  ended. 

Step  6:  If  editing  does  not  occur,  then  the  Kalman  filter  gain  vector, 
Kn+} ,  is  computed: 

*n+l  =  1  - / sr  *  ^n+1  *  ^n  +  1 

The  gain  vector  is  the  measure  of  how  much  the  residual 

r  =  9n+1  -  9n+1 

should  affect  the  state  vector. 

NOTE;  If  a  Kalman  filter  were  designed  to  incorporate  the  three  meas¬ 
urements  (range,  azimuth,  elevation)  from  a  given  laser  at  the  same 
tune,  the  result  would  be  the  same  as  incorporating  them  one  at  a  time; 
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p  . 

but  the  quotient  l./s£  in  the  gain  computation  would  be  a  3x3  matrix 


inversion . 


Step  7:  The  state  vector  estimate  Xn+i  is  now  corrected  optimally  to 
include  the  information  received  from  measurement  6n+i: 


*n+l  =  ^n+l  "  *n+l 


where 


r  =  8n+l  '  8n+l 

Step  8:  The  covariance  matrix  estimate  Pn+j  is  corrected: 

Pn+1  =  ^n+1  "  ^n+1  *  ^n+1  *  ^n+1 

At  this  point,  a  summary  of  the  key  equations  from  the  above  steps 
will  be  written  for  ease  of  reference.  The  equations  assume  tnai:  a 
measurement  has  been  received  at  tune  tn+i  and  that  the  transition  matrix 
to  go  ^t  seconds  from  time  tn  to  tn+i  is  available. 

^n+1  =  ^n»  n+1  *  *n 

^n+1  =  $n»  n+1  *  v^n  *  Q)  *  $n>  n+i 
A  * 

8n+l  =  ^n+1  *  xn+l 
A 

r  =  8  n+1  '  8n+l 

S?  =  Hn+1  *  Pn+1  *  Hn+i  +  [o  (i,  j)}2 

r/sr  =  edit  check 

Kn+1  =  1/S?  *  Hn+1  *  Pn+1 

A 

*n+l  =  ^n+1  "  *Wl  *  r 

^n+1  =  ^n+1  "  ^n+1  *  ^n+l  *  ^n+l 

5.  RAUCH-TUNG-STRIE3EL  SMOOTHER 


The  smoothing  equations  use  a  backwards  recursion  involving  the 
filter  estimates  for  the  state  vector  and  covariance  matrix  to  correct 
earl ier  es timates  based  on  information  received  from  later  measurements. 
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To  accomplish  smoothing,  the  following  information  must  be  recorded 

during  the  Kalman  filtering  of  the  measurements: 
a 

Xn  -  predicted  state  vector  at  time  tn 

X^  -  corrected  state  vector  using  all  measurements  up  to  time  tn 

Pn  -  predicted  covariance  matrix 

P  -  corrected  covariance  matrix  using  all  r  measurements  up  to 

time  tn 

S  ,  ,«  -  transition  matrix  expressing  process  effect  on  X_  to  go 

n  n+i  from  tn  to  tn+1 

The  smoothing  equations  correct  the  state  vector  and  covariance  matrix 
at  time  t  based  on  the  difference  between  their  projection  to  time  tn+1 
and  the  smoothed  values  at  time  tp+1.  Suppose  the  last  time  for  which 
data  was  available  for  the  Kalman  filter  was  tp+1  and  that  all  measure¬ 
ments  have  been  filtered  through  time  tp+j.  Let  Xn+j,  Pp+1  denote 
smoothed  estimates  for  the  state  vector  and  covariance  matrix  at  time  tp+1. 
The  first  step  In  smoothing  is: 

*n+l  '  Xn+1 

Vl  =  Pn+1 

To  obtain  the  smoothed  estimates  at  time  tp>  first  compute  the  gain 
matrix  G: 

15  ■  p„  *  ^  n+l  [Vl]"' 

where 


n 

-T 

^n  n+l 


N" 


-  filtered  estimate  of  covariance  matrix  at  time  t 

n 

-  transpose  of  the  transition  matrix  to  project  the  state 
vector  from  time  tn  to  time  tp+j 

-  inverse  of  the  covariance  matrix  predicted  for  time  tn+j 
based  on  Pn,  the  process  noise  from  tn  to  tn+1  and 

the  transition  matrix  Sn,  n+j 


A 

NOTE:  GPSBET  recomputes  Sn,  n+^  and  Pn+  during  the  smoothing  pass  by 
recording  At  -  tn+i  -  tn  and  the  deweighting  matrix  Q  during  the 
filtering  pass.  This  saves  on  the  amount  of  intermediate 
reading  and  writing  since  the  off  diagonal  elements  of  Q  are  zero 
and  specifying  At  is  all  that  is  required  to  construct  Sn,n+] 

(see  Section  4.3).  The  smoothed  estimate  of  the  state  vector 
at  time  tn  is: 

^n  *  *n  +  l"’  *  (*n+l  '  *n+l) 

The  smoothed  covariance  matrix  at  time  tn  is: 

Pn  =  pn  +  G  *  (Pn+1  -  Pn+l)  *  GT 

This  backwards  recursion  continues  to  the  initialization  time  of 
the  Kalman  filter.  If  a  reinitial ization  occurred  during  filtering 
then  the  smoothing  process  starts  over  just  prior  to  reinitiali¬ 
zation. 

During  the  smoothing  pass,  an  improved  attitude  matrix  is 
computed  from  the  smoothed  velocities  and  acceleration  (unless 
IMU  measurements  were  used). 

Also  during  the  smoothing  pass,  measurement  residuals  are 
computed  and  recorded  on  magnetic  tape.  Measurement  residuals 
are  computed  by  converting  the  smoothed  estimates  of  cartesian 
coordinates  to  the  measurement  domain  as  described  in  4.3.  The 
differences  between  these  predictions  and  measurements  (which  are 
passed  to  smoother  with  the  state  vector  and  covariance  matrix 
estimates)  are  computed. 
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6.  OUTPUT  AVAILABLE 

Output  from  GP5BET  is  in  the  form  of  a  line  printer  listing  and  an 
output  tape.  The  output  tape  has  two  files  for  each  trajectory.  The  first 
file  contains  all  the  program  control  information  such  as  options  selected, 
instrument  locations,  a  priori  information  input  to  the  Kalman  filter, 
etc.  The  second  file  is  a  data  file.  The  line  printer  listing  contains 
the  control  information  in  the  same  form  as  the  first  file  on  the  output 
tape.  The  data  output  on  the  line  printer  is  controlled  by  input  param¬ 
eters  . 

Line  printer  output  during  the  Kalman  filter  pass  can  include  any  or 
all  of  the  following  data  types: 

a.  state  vector  estimates; 

b.  one-sigma  uncertainties  in  state  vector  estimates; 

c.  full  covariance  matrix  estimates; 

d.  measurement  residuals; 

e.  summaries  of  measurement  residuals  over  selected  intervals; 

f.  edit  messages  when  editing  occurs; 

g.  identification  of  retro  switches. 

Data  types  f.  and  g.  are  printed  whenever  they  occur.  Types  a.  through  e. 
are  printed  at  selected  time  intervals  and/or  after  filtering  a  particular 
number  of  measurements.  The  time  intervals  and  points  between  output  are 
selected  by  changing  input  parameters. 

During  the  smoothing  pass  line  printer  output  consists  of  the  state 
vector  and  sigma  vector  (1-sigma  uncertainties  in  state  vector  estimate) 
as  computed  during  the  filter  pass  and  as  computed  in  the  smoothing  pass 
plus  differences  in  the  state  vector  and  the  F-ratio  of  the  variances. 
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A  complete  description  of  the  output  tape  from  GPSBET  is  contained 
in  Appendix  B.  A  sample  of  the  line  printer  output  is  in  Appendix  C. 

7.  VALIDATION 

The  GPSBET  program  has  undergone  extensive  validation  to  ensure  tnat 
the  Kalman  filter  and  smoothing  equations  are  implemented  correctly  and 
that  the  a  priori  information  is  incorporated  properly.  The  validation 
was  accomplished  in  two  stages.  First,  YPu  tested  and  debugged  the  program 
to  the  extent  that  computation  of  -  trajectory  from  a  single  PATS  laser 
compared  favorably  with  YPG's  previous  reduction  programs.  The  next  stage 
of  validation  was  called  the  end-around-check.  The  end-around-check  was 
a  joint  effort  of  The  Aerospace  Corporation  and  YPG.  The  Aerospace  Cor¬ 
poration  was  under  contract  to  the  Space  and  Missile  Systems  Organization  (SAMSO) 
of  the  US  Air  Force.  SAMSO 's  interest  in  YPG's  BET  program  stenmed  from 
the  NAVSTAR  Global  Positioning  System  (GPS)  project.  GPS  intended  to 
(and  does  now)  use  GPSBET  output  as  a  reference  trajectory  for  comparison 
with  the  NAVSTAR  solution  for  position  and  velocity. 

The  end-around-check  included: 

a.  The  Aerospace  Corporation  generated  simulated  trajectories 
similar  to  the  flight  paths  to  be  flown  during  the  GPS  project.  Three 
laser  locations  were  simulated  and  the  X,Y,Z  coordinates  of  the  simulated 
trajectories  were  converted  to  range,  azimuth,  elevation  from  each  of  the 
three  laser  locations.  Random  number  generators  were  used  to  generate 
noise  which  was  added  to  the  measurements.  A  small  percentage  of  the 
measurements  were  changed  to  random  numbers  to  simulate  data  dropouts. 
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b.  The  simulated  laser  measurements  were  transmitted  to  YPG  via 
magnetic  tape  in  the  same  format  as  actual  laser  measurements  were  being 
recorded  at  YPG. 

c.  YPG  personnel  processed  the  simulated  laser  measurements 
with  the  GPSBET  program,  producing  an  output  tape  and  line  printer  listing. 

d.  The  output  tape  from  GPSBET  was  sent  to  The  Aerospace 
Corporation. 

e.  The  Aerospace  Corporation  differenced  YPG  *  s  output  tape 
with  the  true  trajectory  from  which  the  measurements  were  generated. 

f.  The  Aerospace  Corporation  analyzed  the  differences  between 
truth  and  the  estimates  made  by  the  GPSBET  program. 

g.  The  Aerospace  Corporation  produced  reports  of  the  differences, 
their  analysis  of  the  differences  and  suggestions  for  changes  to  GPSBET. 

h.  The  above  cycle  was  repeated  until  The  Aerospace  Corporation, 
SAMSO  and  YPG  were  confident  that  GPSBET  was  filtering  and  smoothing  data 
as  well  as  Kalman  filter  theory  allows  (subject  to  some  time  and  money 
constraints) . 

i.  All  of  the  above  steps  were  repeated  with  the  same  trajec¬ 
tories  but  using  simulated  inertial  platform  measurements  in  addition  to 
laser  measurements. 

The  results  of  the  end-around-checks  are  summarized  in  Tables  1  and 
2  which  contain  the  RMS  of  the  difference  between  truth  and  GPSBET  esti¬ 
mates  for  two  trajectories  called  I  and  II  and  for  two  types  of  data  A 
and  B.  Trajectory  I  was  a  gentle  .5g  turn,  whereas  trajectory  II  was  a 
figure-8  with  a  4.5g  turn  on  one  end  and  a  3 . 5g  turn  on  the  other  in 
addition  to  Ig  vertical  climbs  and  descents.  Cases  IA  and  IIA  were  using 


only  simulated,  laser  data  whereas  Cases  IB  and  1 1 B  were  using  simulated 
data  from  3  lasers  plus  simulated  inertial  platform  data.  Table  1  summa¬ 
rizes  how  well  the  filter  only  estimates  compared  with  truth,  whereas 
Table  2  summarizes  tne  results  of  filtering  and  smoothing.  Timing  results 
in  Section  8.3  show  how  expensive  the  smoothing  pass  is. 


TABLE  1* 


FILTER  ONLY 

RMS  OF  DIFFERENCES  WITH  TRUTH 


IA 

IB 

I  IA 

I  IB 

UNITS 

X 

m 

.07 

.02 

.40 

.24 

Y 

m 

.03 

.01 

.30 

.17 

Z 

m 

.03 

.01 

.48 

.24 

X 

m/sec 

.18 

.01 

1.40 

.45 

Y 

m/s  ec 

.08 

.01 

1.01 

.33 

Z 

m/sec 

.02 

.00 

1.23 

.28 

x’ 

m/sec2 

.27 

.07 

3.24 

1.76 

y‘ 

m/sec2 

.11 

.03 

2.13 

1.22 

z 

m/sec2 

.02 

.00 

1.99 

.91 

*  These  are  results  of  analysis  by  The  Aerospace  Corporation 
and  were  taken  from  DPDWG  Technical  Memorandums  prepared  by 
H.  Bernstein,  The  Aerospace  Corporation. 


TABLE  2* 


FILTER  AND  SMOOTH 
RMS  OF  DIFFERENCES  WITH  TRUTH 


IA 

IB 

I  IA 

lie 

UNITS 

X 

m 

.02 

.01 

.23 

.ii 

Y 

m 

.01 

.01 

.17 

.08 

Z 

m 

.01 

.01 

.28 

.11 

X 

m/sec 

.01 

.00 

.38 

.10 

Y 

m/sec 

.01 

.01 

.30 

.09 

Z 

m/sec 

.01 

.00 

.40 

.09 

x' 

m/sec2 

.03 

.02 

1.13 

.64 

Y* 

m/sec2 

.02 

.01 

.89 

.50 

Z 

m/sec2 

.01 

.00 

.90 

.43 

*  These 

are  results 

of  analysis  by  The  Aerospace  Corporation 

and  were  taken  from 

OPDWG 

Technical 

Memorandums 

prepared  by 

H.  Bernstein,  The  Aerospace  Corporation. 
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8.  PERFORMANCE  USING  ACTUAL  DATA 


The  performance  of  GPSBET  using  actual  data  is  monitored  very  closely 
since  VPG  customers  depend  on  the  absolute  accuracy  of  the  reference 
trajectory  that  GPSBET  produces.  YPG  maintains  extensive  histories  of 
laser  calibration  coefficients,  monitors  measurement  noise,  computes 
process  noise,  studies  measurement  redisuals,  plots  inter-laser  differ¬ 
ences,  plots  single  laser  solutions  against  GPSBET  solutions,  etc.,  to 
ensure  that  GPSBET  output  is  as  good  as  the  measurements  allow.  In 
addition  to  the  internal  checks  mentioned  above,  YPG  invites  comparisons 
with  instrumentation  and  reduction  procedures  external  to  the  laser 
systems  and  YPG  software.  Sections  8.1  and  8.2  describe  some  examples 
of  the  internal  and  external  comparisons  and  analysis  that  YPG  has  done 
and  will  continue  to  do  to  verify  the  accuracy  of  GPSBET  estimates. 

8.1  INTERNAL  ANALYSIS 

Extensive  analysis  of  GPSBET  accuracy  must  be  performed  since  range 
users  expect  0.5  meter  accuracy  in  position  estimates.  The  most  important 
requirement  for  attaining  this  accuracy  Is  pre-  and  post-mission  calibra¬ 
tion  of  all  lasers  to  update  calibration  histories  and  allow  analysis  of 
potential  problems.  Another  requirement  is  complete  meteorological  data 
for  accurate  atmospheric  corrections.  A  priori  information  about  measure¬ 
ment  noise  and  process  noise  must  be  verified.  YPG's  real-time  data 
collection,  computation  and  display  facility  provides  a  means  of  acquiring 
the  information  required  for  an  accurate  trajectory  estimate. 

After  a  real-time  mission,  the  standard  procedure  is  to  generate 
quick-look  plots  of  the  trajectory  data.  These  plots  include  at  least: 
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X  vs  Y 

Y  vs  7. 

X,Y,Z  vs  time 

X,Y,Z  vs  tin* 

RTE  source  vs  time 

The  qui..c -look  plots  are  generally  available  within  4  hours  after  a  2-hGur 
mission.  These  plots  can  be  used  to  adjust  the  a  priori  information  for 
GPS8ET. 

In  many  cases,  the  ReaV-Tiine  Estimate  (RTE)  of  the  trajectory  of  or. 
aircraft  is  more  than  adequate  to  meet  customer  accuracy  requirements. 

In  particular,  when  a  customer  uses  only  a  single  laser  for  data  acquisi¬ 
tion,  then  the  RTE  is  almost  as  accurate  as  the  GPSBET  for  position  and 
velocity.  The  Kalman  filter  and  smoother  output  is  however,  less  noisy 
and  contains  good  acceleration  estimates.  Figure  1  shows  an  X  vs.  Y  anc 
Y  vs  Z  plot  of  a  segment  of  a  helicopter  trajectory  for  which  some  special 
analysis  was  done.  A  single  laser  (Site  7  laser)  was  forced  to  be  used 
for  the  RTE  (the  source  selection  algorithm  used  in  real  time  was  manually 
overriden)  and  the  GPSBET  was  run  using  data  only  from  Site  7  laser. 
Figures  2,  3,  4  show  the  X,  Y,  Z  differences  between  the  RTE  and  BET 
processing  when  both  used  data  from  Site  7  laser  only.  Differences  are 
attributable  to  laser  data  dropouts,  noise  in  the  estimates,  some  differ¬ 
ences  in  the  level  of  sophistication  in  atmospheric  correction  and  compu¬ 
tational  errors  due  to  roundoff  and  truncation. 

For  a  200-second  segment  of  the  trajectory  shown  in  Figure  1,  a 
three  laser  BET  was  computed.  The  three  laser  BET  was  differenced  with 
each  of  the  individual  laser  solutions  as  computed  in  real  time.  Plots 
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of  the  differences  appear  in  Figures  5-13.  Figures  5,  6,  7  are  the  X,  Y,  Z 
differences  between  the  Site  7  laser  solution  computed  in  real  time  and 
the  three  laser  BET.  Figures  8,  9,  10  are  the  X,  Y,  Z  differences  between 
the  Site  9  laser  solution  computed  in  real  time  and  the  three  laser  BET. 
Figures  11,  12,  13  are  the  X,  Y,  Z  differences  between  the  Site  12  laser 
solution  computed  in  real  time  and  the  three  laser  BET. 

Contained  on  the  output  tape  from  the  three  laser  BET  run  made  on 
the  200-second  segment  of  trajectory  referred  to  above  were  the  measure¬ 
ment  residuals  for  the  three  lasers.  Measurement  residuals  are  computed 
by  transforming  the  best  estimate  of  X,  Y,  Z  at  each  measurement  time  to 
range,  azimuth,  elevation  from  each  laser  and  differencing  them  with  the 
measurements.  Figures  14-22  contain  plots  of  the  range,  azimuth,  eleva¬ 
tion  residuals  for  Sites  7,  9,  12.  Residuals  plots  show  what  the  magni¬ 
tude  of  the  measurement  noise  is  and  nelp  to  identify  calibration  problems. 
Notice  that  noise  on  the  range  measurements  is  about  .5  meter  and  on 
azimuth  and  elevation  is  about  .1  milliradian.  The  residuals  plots  in 
Figures  14-22  were  generated  to  analyze  a  suspected  accuracy  problem  and 
provided  additional  verification  of  a  Site  12  azimuth  encoder  problem 
which  was  subsequently  corrected. 

8.2  COMPARISON  WITH  OTHER  INSTRUMENTATION  AND  PROCESSING 

GPSBET  trajectory  estimates  have  been  compared  with  many  other  types 
of  processing  using  the  laser  instrumentation  and  other  YPG  instrumenta¬ 
tion  such  as  cinetheodol i tes ,  PLS  (Position  Location  System,  which  is  an 
RMS  II  System  manufactured  by  General  Dynamics)  and  several  radars.  Als  , 
comparisons  were  made  with  an  airborne  navigation  solution  generated  by 
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CIRIS.  The  GPS  project  has  repeatedly  demonstrated  agreement  within 
3.  meters  of  the  GPSBET  solutions.  In  support  of  the  GPS  project,  The 
Aerospace  Corporation  processes  laser  and  INS  data  independently  and  con¬ 
sistently  demonstrates  1  sigma  differences  within  .25  meter  of  YPG's 
GPSBET  trajectory  estimates. 

The  most  conclusive  proof  of  .5  meter  for  the  absolute  accuracy 
(1  sigma)  of  GPSBET  position  estimates  was  derived  from  a  Ballistic 
Camera  Comparison  Test.  The  test  was  conducted  in  January  1977.  USATECOM 
funded  White  Sands  Missile  Range  to  set  up  four  ballistic  cameras  on  YPG's 
Cibola  Range. 

Data  was  collected  for  many  passes  of  an  F-4  jet  through  a  volume  of 
airspace  providing  good  geometry  for  the  ballistic  cameras.  Because  of 
inclement  weather  and  other  data  collection  problems  (a  disadvantage  of 
ballistic  camera  data  collection),  only  eight  passes  of  ballistic  camera 
data  could  be  reduced.  Physical  Sciences  Laboratory  (PSL)  computed  tra¬ 
jectory  estimates  using  the  ballistic  camera  data  and  provided  their  final 
output  to  YPG  on  magnetic  tape.  YPG  processed  tne  data  from  the  three 
lasers  through  the  GPSBET  program.  The  GPSBET  program  made  lever  arm 
corrections  so  that  the  solution  from  GPSBET  was  for  the  position  of  the 
strobe  light  used  in  the  ballistic  camera  reduction.  Figures  23-26  are 
X  vs.Y  plots  showing  the  ground  track  of  the  first  four  trajectories 
analyzed.  The  F-4  target  aircraft  was  flying  at  a  15000  foot  altitude 
for  all  passes. 

The  ballistic  camera  solution  had  a  1  sigma  uncertainty  of  about 
.2  meter  near  the  center  of  each  pass.  Near  the  ends,  the  solution  de¬ 
graded  to  a  two-camera  solution  for  some  passes  and  the  uncertainty 
increased  to  more  than  1.0  meter. 
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Differences  between  the  GPSBET  solutions  using  three  lasers  and  the 
ballistic  camera  solutions  were  about  .5  meter.  Table  3  contains  the 
rms  (root  mean  square)  of  the  differences  in  X,  Y,  Z  components  for  each 
of  the  eight  passes  for  which  ballistic  camera  data  was  available. 

Figures  27-38  are  plots  of  the  point-by-point  differences  between  GPSBET 
and  ballistic  camera  solutions  for  X,  Y,  Z  components  in  passes  1-4. 

8.3  COMPUTER  TIME 

The  GPSBET  program  is  Implemented  on  an  IBM  7094/7044  DCS  computer. 
Timings  were  made  using  various  types  of  instrumentation.  Table  4 
contains  the  results  of  the  timings  along  with  the  approximate  uncertain¬ 
ties  estimated  for  the  position  and  velocity  components.  The  run  time 
column  represents  multiples  of  the  length  of  the  trajectory  being  esti¬ 
mated.  These  times  explain  why  the  GPSBET  program  is  not  usually  run 
over  an  entire  mission  but  over  smaller  time  segments  determined  from 
quick-look  plots  using  RTE  data. 


TABLE  3. 


RMS  OF  DIFFERENCES  IN  X.  V.  Z 
LASER  CALIBRATIONS  ADJUSTED  FOR  TEMPERATURE  INVERSION 

(«et«rs) 

Pass  1  Pass  Z  Pass  3  Pass  4  Pass  5  Pass  6  Pass  7  Pass  2 


00.00  SaiO.QO  5&20.00  51330.00  5340.00 
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9.  SUMMARY 


USAYPG  acquired  three  laser  tracking  instruments  in  1973-1974. 
Accompanying  this  acquisition  was  the  implicit  assumption  that  USAYPG 
personnel  would  develop  standard  operating  procedures,  design  hardware 
and  design  and  develop  software  to  take  full  advantage  of  these  new 
tracking  systems.  Since  it  turned  out  that  the  lasers  had  much  greater 
capabilities  for  accurate  and  timely  data  acquisition  and  reduction 
than  any  previous  range  instrumentation,  this  implicit  assumption  meant 
devising  calibration  procedures  and  trajectory  estimation  proceoures  at 
or  beyond  the  state-of-the-art.  From  wnence  came  the  requirement  for 
the  GPSBET  program. 

As  described  in  this  paper,  the  GPSBET  program  computes  optimal 
(or  as  nearly  optimal  as  the  truth  of  the  assumptions  allow)  estimates 
of  position,  velocity  and  acceleration  utilizing  range,  azimuth  and 
elevation  measurements  from  up  to  tnree  laser  trackers.  Also,  GPSBET 
accepts  inertial  platform  measurements  in  the  form  of  sensed  velocity 
components  and  gimbal  angles. 

The  GPSBET  program  accepts  control  information  from  cards  and  data 
from  magnetic  tape.  It  processes  measurements  through  Kalman  filter 
equations  followed  by  an  optional  smoothing  pass. 

GPSBET  is  implemented  on  an  IBM  7094/7044  Direct  Couple  System. 

On  that  computer,  the  Kalman  filtering  of  data  from  three  lasers  takes 
about  four  times  real  time  (i.e.,  a  15-minute  segment  of  data  takes  1 
hour  of  computer  time).  Filtering  plus  smoothing  takes  about  nine 
times  real  time. 

The  accuracy  of  position  estimates  made  by  GPSBET  using  data 
from  three  lasers  has  been  demonstrated  to  be  about  .5  meter 
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Appendix  A 


BET  INPUT  CONTROL  PARAMETERS 

This  appendix  lists  and  describes  all  parameters  that  can  be 
input  to  the  BET  program  from  cards  as  of  15  April  1976. 


r'ard i  iQ ter  NdH.e 

NJ0B  -  -  -  - 

NT RAJ  -  -  -  - 

GC0RSQ  -  -  -  - 

GC0RS1  -  -  -  - 

GC0RS2  -  -  -  - 

GC0RS3 
RETLEV 

IMULEV  -  -  -  - 

GPSLEV 

XS  -  -  -  - 

NI 

INSTYP 


INUNIT  -  -  -  - 
I0BET 

IMU  -  -  -  - 
ISMUTH  -  -  -  - 
ILEV  -  -  -  - 


Descr ption 
number  of  jobs  to  do 
number  of  tra jectories  per  job 

geodetic  coordinates  of  origin  of  reference  coordinate 
system 

geodetic  coordinates  of  laser  1 
geodetic  coordinates  of  laser  2 
geodetic  coordinates  of  laser  3 
body  coordinates  of  retroref lector 
body  coordinates  of  IMU 
body  coordinates  of  GPS  antenna 
state  vector 
number  of  Instruments 

type  of  each  Instrument 

1  =  laser 

2  ■  cine 

3  *  range  only 

4  ■  sensed  velocities 

5  *  glmbal  angles 

logical  unit  number  to  read  COT 
output  unit  for  BET  output  tape 
option  to  use  IMU  data 
option  to  make  smoothing  pass 
option  to  make  lever  arm  corrections 
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<r. 


Por<»M-'ter  ,\ame 


Description 


-  -  -  -  option  to  rewind  CDT  after  a  run 


IRDE0F  -  -  - 


INSTRU  -  - 


j  kA  T  t 


IXSb  • 


MAXPTS  - 


I  FI LES  - 


option  to  read  to  an  E0F  on  CDT  after  filtering 
a  given  section  of  data 

Instrument  to  De  forced  as  automatic  Initialization 
or  reinitialization  of  the  BET 

array  controlling  nodules  for  output  of  various 
types  of  data  to  line  printer  or  tape 

max  points  Detween  stute  vector  output  during 
f i 1 teriny 

max  points  between  sigma  vector  output  during 
f i 1 tering 

max  points  between  printout  and  restart  of  rms  and 
mean  of  residuals  report  during  filtering 

max  points  between  printout  of  actual  residuals  for 
a  single  time  from  each  instrument 

max  points  between  printout  of  covariance  matrix 
(actually  has  sigmas  on  diagonal  and  correlation 
coefficients  in  off  diagonal  elements) 

maximum  number  of  observations  to  filter  on  this 
trajectory 

number  of  files  to  skip  on  the  CDT  prior  to  reading 
the  data  for  this  trajectory 


I  REF  ---  -  option  on  refraction  corrections 

XlAMDA  -  -  -  -  wavelength  of  lasers 

HSCL  -  -  -  -  scale  height  parameter  used  in  current  refraction 

mo  Gel 

NIMU  -  -  -  -  number  of  the  receiver  to  use  sensed  velocities 
from  (number  refers  to  location  on  the  CDT) 

NGMBL  -  -  -  -  number  of  the  receiver  to  take  gimbal  angles  from 

M  ....  array  containing  options  for  prefiltering  or  skipping 

0  =  skip 

1  =  constant  fit 

2  *  linear  fit 

3  =  quadratic  fit 

N  -  -  -  -  number  of  points  to  use  in  prefilter 


108 


Description 


Parameter  Name 

PLDNSD  -  -  -  -  number  of  multiples  of  residual  standard  deviation 
to  use  while  editing  prior  to  the  pre-filter 

RAWDT-  -  -  -  -  time  increment  between  raw  data  for  each  instrument 

FILDT  -  -  -  -  time  Increment  desired  between  measurements  going 

through  the  Kalman  filter 

EDIT  -  -  -  -  number  of  standard  deviations  of  residuals  to  edit 
on  during  the  Kalman  filtering  (after  pre-filtering) 

IEDLMT  -  -  -  -  number  of  edit  messages  to  print  out  during  this  run 

SD  -  -  -  -  standard  deviations  of  noise  on  measurements  of  each 

tape  and  from  each  instrument 

QI  -  -  -  -  standard  deviation  used  for  deweighting  the  covariance 

matri x 

PI  -  -  -  -  initial  values  for  square  root  of  the  diagonal  of 
the  covariance  matrix  (i.e.>  initial  uncertainty  in 
state  vector) 

TEMP  -  -  -  -  temperature  in  degrees  Kelvin  at  each  laser  site 

PRES  -  -  -  -  pressure  in  millibars  at  each  laser  site 

RUN  I D  -  -  -  -  72  character  run  identification 

TEST  I  -  -  -  -  20  character  test  item  identification 

IW0N  -  -  -  -  12  character  work  order  number  or  other  identification 

of  project  to  be  charged 

PRjJJ  -  -  -  -  12  characters,  for  project  engineer's  name 

DATE  -  -  -  -  12  characters  for  date  of  test 

RDATE  -  -  -  -  12  characters  for  date  of  this  report 

NHI 

NMI 

NS  I 

MSI  -  -  -  -  start  and  stop  times  for  using  measurements  from 

NHF  |  -  -  -  -  various  instruments  in  the  filter 

NMF 

NSF 

MSF 
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1 


Parameter  Name 


description 


NHIESLT 
NK1BLT 
NSIBCT 
MS  IBtT 
NHFBLT 
NMFBtT 
NSFBFT 
MSFBtT 


start  and  stop  times  for  filtering  on  this 
trajectory 
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Appendix  B 

BET  OUTPUT  TAPE 

For  each  trajectory  filtered,  there  will  be  a  BCD  file  containing 
control  information  and  a  second  file  containing  results.  This 
appendix  includes  a  description  of  the  second  file.  The  first  file  for 
each  trajectory  contains  a  copy  of  the  control  report,  a  sample  of 
which  is  in  Appendix  C  of  this  paper. 


YUMA  PROVING  GROUND 

DATA  PROCESSING  DEVELOPING  WORKING  GROUP 


DPDWG  REPORT  NO.  2 
BET  OUTPUT  RECORD  FORMAT 


18  December  1975 


ABSTRACT 

The  agreed-upon  record  format  in  which  the  output 
from  the  Best  Estimate  of  the  Trajectory  (BET)  program 
will  be  written  is  presented;  the  information  to  be 
contained  in  the  BET  output  is  defined  by  this  format. 

Changes  to  this  format  must  be  documented  in  sub¬ 
sequent  DPDWG  reports. 
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BFT  OUTl'l'T  RECORD  FORMAT 


The  output  from  the  B F T  (Best  Estimate  of  the  Trajectory) 
program  will  be  written  by  the  IBM  7014/7094  DCS  in  JG-bit  word* 
on  7-track  magnetic  tape.  The  output  from  each  run  will  consist  of 
two  fijos: 

File  1.  The  first  file  will  be  written  in  BCD,  containing 
all  of  the  information,  and  in  identical  form,  of 
the  int roductory  portion  ("control  report")  of  the 
BET  printout.  The  first  file  will  thus  contain 
date  and  time  of  the  test,  identification  and  location 
of  range  sensors,  and  all  input  parameters. 

File  2.  The  second  file  will  contain  the  results.  It  will  be 
written  in  binary  words,  with  one  record  each 
0.  2  seconds,  in  order  of  increasing  time,  and  one 
physical  record  per  logical  record. 

If  a  physical  BET  output  tape  contains  the  output  from  several 
runs,  there  will  be  two  files--the  control  report  and  the  results  file-- 
for  each  run. 

Table  I  summarizes  the  record  format  which  will  be  employed  in 
the  second  (results)  file. 

Rcpr<  '.entation  of  Words.  A  word  in  the  results  record  will  be  written 
by  the  7044/7094  DCS  as  a  result  of  FORTRAN  (binary)  write  commands, 
either  as  an  integer  word  or  as  a  floating  point  word. 

Fixed  Point  ( Integ  r  r  )  Word 


TuTd 


The  sign  of  the  number  is  contained  in  bit  36.  A  "0‘  signifies  a 
positive  number  and  ,i  ’1"  signifies  a  negative  number.  The 
magnitude  of  the  number  is  in  bits  35  (most  significant)  through  1 
(least  significant),  with  the  binary  point  positioned  to  the  right  of 
bit  1 . 


Floating  P  urn  M  ‘-uncle  Precision  Real)  Wo  rd 
. - -.CJurav  lcxi  _ _  Fraction 


36  35  28  27  1 

i  he  sign  of  the  number  is  contained  in  bit  36.  In  bits  35  through 
28  i<>  Uk'  binary  exponent  mcreiscd  by  128  (decimal)  200  (octal) 
10000000  (binary).  In  bits  27  through  1  is  the  fraction  (mantissa), 
with  the  hiniry  point  positioned  between  bits  28  and  27. 
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Aircraft  attitude  matrix  N/A  Real 


FF.T  OUTPUT  RECORD  FORMAT  {CONT'D> 


i  m.i 


(  \V  urn  II 


Wcn'il  1  givi  s  the  turn-,  in  seconds,  of  L’T  (''Universal  time, 
or  (,M1  ).  Since  local  time  at  Yuma  Proving  Ground  is  Mountain  Standard 
1  line  (MSI  ),  tlie  tune  in  word  1  will  be  greater  by  7  hr  =  25200  sec  than 
local  time. 

Source  (Word  2) 

Word  2  indicates: 

(1)  Whether  the  BET  is  the  result  of  a  filter  pass  only,  or  the 
result  of  a  smoother  pass. 

(2)  Which  sensors  contributed  to  the  solution  in  the  filter  (noi 
smoother)  pass.  A  sensor  will  be  indicated  as  contributing 
if  it  had  meaningful  data  at  the  time  its  last  measurement 
set  was  due  and  at  least  two  of  three  of  the  measurements 
were  accepted  (not  edited  in  the  filter). 

Table  11  gives  the  various  values  which  word  2  (source)  Can  take  on. 

BFT  'state  Vector  (Words  3-26) 

The  BET,  based  upon  measurements  from  one  to  three  laser 
trackers  and,  on  some  flights,  on  inertial  measurement  unit,  is  gen¬ 
erated  by  a  Kalman  filter  with  an  optional  smoother.  The  program 
develops  estimates  of  the  host  vehicle  position,  velocity,  and  acceleration, 
of  the  range,  azimuth,  and  elevation  biases  of  each  laser  tracker,  and  of 
the  orientation  and  velocity  biases  from  the  inertial  navigation  system  (INS). 
These  quantities  are  referred  to  as  the  "state  vector,"  which  will  contain 
from  12  to  24  elements,  depending  upon  the  number  of  sensors  employed. 
The  program  also  computes  the  covariance  matrix  containing  the  expected 
errors  of  the  state  vector  parameters  and  the  expected  correlations 
among  the  errors. 

Normally,  the  BET  will  be  (lie  result  of  a  smoother  pass,  and  the 
state  vector  (words  3-2u)  will  be  generated  by  the-  smoother  pass.  Should 

It  is  anticipated  that  Y  PG  range  time  will  he  expressed  in  UT  for  GPS 
testing.  Should  raiute  time  still  be  expressed  in  MST,  the  correction 
will  be  performed  lay  software. 


116 


tamle  ir 


souuc: i  identification  ron  ret  (format  word  Z) 


SOURCE 

VALUE 

Filter  Pass  Smoother  Pass 

GPS  (Telemetered  Data) 

-999999 

999999 

No  data  (extrapolated) 

-99999 

99999 

Laser  1 

-1 

1 

Laser  2 

-2 

2 

Laser  3 

-3 

3 

PLS 

-4 

4 

Lasers  1  &  2 

-12 

12 

Lasers  1  k  3 

-13 

13 

Lasers  2  &  3 

-23 

23 

Lasers  1  &  2  &  3 

-123 

123 

IMU 

-1000 

1000 

IMU  +  Laser  1 

-1001 

1001 

IMU  +  Laser  2 

-1002 

1002 

IMU  4  Laser  3 

-1003 

1003 

IMU  +  PLS 

-1004 

1004 

IMU  +  Lasers  1  k  2 

-  1 0 1  2 

1012 

IMU  3  Lasers  1  k  3 

-1013 

1013 

IMU  +  Lasers  2  &  3 

-1023 

1023 

IMU  -f  Lasers  1  Kt  2  fci  3 

-1123 

1123 
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t'.t  I'l.J  be  generated  by  the  lilu  r  pass  only,  the  state  vector,  of  course, 
will  i  •.line  from  th.it  pass. 

T  r.i  |i  t  t  or  y 

K  t-fe  ren<  e  1  *01111 .  Words  3-11  contain  the  position,  velocity,  and 
acceleration  of  the  point  for  which  the  BET  is  calculated.  If  the  BET  is 
bused  entirely  upon  laser  data,  the  calculation  point  will  be  the  laser 
ret  rorcflector.  If  the  BET  is  ba sod  upon  INS/1MU  plus  laser  data  or 
upon  INS/IMU  plus  PLS  data,  the  calculation  point  will  be  the  gimbal 
center  of  the  IMU. 

Coordinate  System.  The  trajectory  will  be  referred  to  a  cartesian 
coordinate  system  ("IRCC  coordinate  system")  which  will  be  centered  at 
the  IRCC  with  the  x,  y,  z  axis  directed  in  the  local  east,  north,  vertical 
directions,  respectively. 

Position.  Words  3,  4,  and  5  contain  x,  y,  and  z,  the  coordinates 
of  the  reference  point  (  retroreflector  or  IMU);  position  coordinates  will 
be  expressed  in  meters. 

Velocity.  Words  6.  7,  and  8  contain  x,  y,  and  i,  the  velocity 
components  of  the  reference  point,  expressed  in  meters /sec. 

Acceleration.  Words  9,  10,  11  contain  x,  y,  z.  the  acceleration 

2 

components  of  the  reference  point,  expressed  in  meter/sec  . 

Laser  Bins  States 

kascr  1  Word  12  contains  the  filter  estimate  of  the  range  bias  of 
Laser  1,  expressed  in  meters. 

Word  13  contains  the  filter  estimate  of  the  azimuth 
bias  of  Laser,  expressed  in  radians. 

Word  14  contains  the  filter  estimate  of  the  elevation 
bias  of  Laser  1,  expressed  in  radians. 

Laser  2  Words  19-17  <ontainthe  filter  estimates  of  the  range, 
azimuth,  and  elevation  biases  of  Laser  2,  expressed, 
respectively,  in  meters,  radians,  and  radians.  Should 
only  one  laser  have  been  employed,  words  15-17  will 
contain  zero. 
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Laser  3 


Words  1  3-20  contain  the  filter  estimates  of  the  range, 
azimuth  and  elevation  biases  of  Laser  3,  expressed, 
respectively,  in  meters,  radians,  and  radians.  Should 
only  one  or  two  lasers  have  been  employed,  words 
lb-20  will  contain  zero. 

INS  Bias  States 

Words  21-26  contain  the  filter  estimates  of  the  INS  bias  state. 

Should  no  IMU  have  been  employed,  words  21-26  will  contain  zero. 

Words  21-23  contain  the  estimated  INS  orientation  biases  about  the 
INS  cast,  north,  and  vertical  axes,  expressed  radians. 

Words  24-26  contain  the  estimated  INS  velocity  biases  in  the  INS 
east,  north,  and  vertical  directions,  expressed  in  meter/sec. 

Standard  Deviations  of  Slate  Vector  Elements  (Words  27-50) 

Words  27-50  contain  the  uncertainty  (standard  deviations)  in  the 
estimated  values  of  the  state  variables,  computed  by  taking  the  square 
roots  of  the  diagonal  elements  of  the  covariance  matrix.  The  uncertainties 
are  in  the  same  order  as  the  state  vector  elements,  and  expressed  in  the 
same  units.  Thus  word  27  contains  the  standard  deviation  of  the  x  position 
coordinate  (word  3)  and  is  expressed  in  meters;  similarly,  word  41  contains 
the  standard  deviation  of  the  elevation  bias  of  Laser  2  (word  17)  and  is 
expressed  in  radians. 

When  a  bias  siatc  word  contains  zero  because  the  sensor  was  not 
employed,  the  corresponding  standard  deviation  will  likewise  contain  zero. 

Aircraft  Attitude  (Words  51-5°) 

In  order  to  compare  the  GPS  solution  with  the  BET,  it  is  necessary 
to  transform  the  BET  position  to  the  position  of  the  GPS  reference  location 
on  the  aircraft.  This  transformation  is 

X  =  X*  +  B  6X, 
r  <  b 

in  winch  B  is  a  tr, information  (rotation)  matrix  and 
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r~rr 


x’ 

r 


position  oi  (<J'v>  r<  firrm  i  luxation,  iri  IKCC  cooruinate s 


X* 

c 


puMiKitj  ui  point  for  whit  h  HKT  was  i  air ulatod ,  in 
IRCC  coordinates. 


relative  location  of  to  PS  reference-  location  to  BET 
Calculation  point,  in  aircraft  body  coordinates. 


The  3x3  matrix  lUl.Jl  is  contained  in  words  51  through  59. 
first  and  second  indices  tlenote  row  and  column,  respectively. 


Word 

Matrix  Element 

51 

B( l .  1) 

52 

B(2.  1) 

53 

B(3.  1) 

54 

B(  1 . 2) 

55 

B(2,  2) 

56 

B(3,  2) 

57 

B(l.  3) 

58 

B(2,  3) 

59 

B(3.  3) 

Covariance  Matrix  (Words  60-65) 

In  words  60-65  are  recorde-d,  in  packed  form,  the  covariance 
elements  corresponding  to  position  states.  The  matrix  is  written  out 
by  half  columns 


,  Word 

Covariance  Element 

60 

P<  1 .  1 ) 

= 

a2 

X 

6  i 

P(  1  .  2)  - 

P(2.1)  = 

^2°x°y 

62 

P(2.  D 

- 

o2 

y 

63 

PO  .  3)  = 

P<3.  1)  = 

P , ,  o  O 
13  x  z 

64 

II 

CO 

OJ 

rT' 

Pt  3 , 2 )  - 

p,,  o  o 

2  3  y  z 

65 

P(  3.  3) 

s 

2 

a 

z 


with  the  p  bome  correlation  coefficn  nts. 

I  ,J 
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Words  66-80.  i  oat. i no ng  the  ritim  atid  values  of  the  state  variables 
wlm  h  were  computed  during  the  filter  pass,  rorrrspond  in  order,  units, 
and  representation  to  words  }-/<•.  Should  the-  BET  be  based  upon  the 
filter  pass  aloiu  ,  words  o6-89  will  all  contain  zero. 

St  a  nd  a  rd  Deviation  of  State  Vector  Elements  from  Filter  Solution  (Words  90-113) 

Words  90-  1  1  3.  containing  the  standard  deviations  of  the  values  of 
the  state  variables  estimated  during  the  filter  pass,  correspond  in  order, 
units,  and  repre sentat ion  to  words  27-50.  Should  the  BET  be  based  upon 
the  filter  pass  alone,  words  90-113  will  all  contain  zero. 

Table  III  summarizes  the  location  in  the  BET  output  record  of  the 
values  of  the  state  vector  elements  and  their  uncertainties,  depending 
upon  whether  the  parameter  was  determined  in  the  smoother  pass  or  in 
the  filter  pass. 

Measurement  Residuals  (Words  ll4-12rO 

The  capability  to  calculate,  by  means  of  a  measurement  model, 
the  value  which  a  measurement  would  give  in  the  absence  of  noise,  is 
inherent  in  a  Kalman  filter.  (The  measurement  model  calculation  is 
based  upon  the  value  of  the  state  vector  elements.  )  The  measurement 
residuals  contained  m  words  114-125  are  the  result  of  applying  the 
measurement  model  to  the  BET  state  vector  and  forming  the  difference: 

(Residual)  ("Measurement")  -  (Measurement  Model) 

For  the  laser  measurements,  the  "measurement"  is  the  value 
produced  by  pre-editing  laser  data  and  pre-filtering  from  20  samples 
per  second  to  5  samples  par  second.  A  separate  laser  residual  is 
calculated  for  each  MET  output  record.  If  a  laser  was  net  employed, 
its  measurement  residuals  will  contain  zero. 

The  data  interval  for  IMl'/INS  data  is  expected  to  be  variable, 
wit  tun  a  pass,  but  to  be  great  it  than  0.2  see.  Consequently,  the  IN'S 
measurement  residuals  will  be  mi  1  tided  if  an  INS  measurement  set 
oe  i  n  r  red  within  the  previous  0,2  sec  i  rid  will  be  zero  otherwise.  If 
an  JMII  was  not  employed,  all  I  NS  nun  su  rement  residuals  will  be  o  ru. 
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LOCATION  OF  STATE  VECTOR  ELEMENTS  AND  UNCERTAINTIES 


APPENDIX  C 


This  appendix  contains  samples  of  the  types  of  data  available  as 
line  printer  output  from  the  GPSBET  program.  The  first  four  pages  are 
an  example  of  the  Kalman  filter  control  report.  This  report  lists  input 
control  parameters  and  a  priori  information  for  the  Kalman  filter.  Next 
are  the  first  three  pages  of  a  typical  run  of  the  Kalman  filter  accepting 
data  from  three  lasers  and  INS  (24  components  in  the  state  vector).  The 
last  two  pages  of  this  appendix  are  copies  of  the  first  two  pages  of  the 
Kalman  smoothing  report.  Note  that  output  from  Kalman  smoothing  is  in 
reverse  chronological  order  (time  backs  up)  since  smoothing  is  accomplished 
by  a  backwards  recursion  on  the  filtered  output. 
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Appendix  D 


This  appendix  contains  plots  of  differences  between  GPSBET 
estimates  for  a  trajectory  and  the  true  trajectory  from  which  simulated 
measurements  were  generated. 
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INTRODUCTION 


1  .0 

The  TAPP*  (Trajectory  Analysis  and  Prediction  Program)  is  designed  as  a 
computational  tool  to  simulate  trajectory  motion  and  tracking  operations,  to 
reconstruct  model  characteristics  responsible  for  observed  behavior  and  to 
perform  error  analyses  by  propagating  uncertainties  in  initial  conditions, 
measurements  and  sensor  models  into  trajectory  position,  velocity  and 
acceleration  at  subsequent  times. 

1.1  Background 

In  mid-1966  the  TRACE-D  (IBM  7094)  orbit  determination  program  was  obtained 
at  AFETR  for  evaluation.  Through  the  remainder  of  1966  and  most  of  1967 
close  contact  was  established  with  the  Aerospace  Corporation  for  improving 
and  using  this  version  of  TRACE. 

Late  in  1967  TRACE  66**  was  made  available  for  use  on  the  CDC  6600  and  3600 
computers.  This  program  was  considerably  improved  over  1RACE-D  and  included 
years  of  experience  gained  from  its  predecessor.  It  was  quoted  by  Aerospace 
to  represent  a  $6M  effort  to  that  point. 

During  the  summer  of  1968  the  TRACE  66  program  was  converted  to  run  on  the 
IBM  360/65  by  RCA  with  several  extensions  aimed  at  improving  its  endoat- 
mosphere  trajectory  solution  capabilities.  This  version  was  renamed  FIAT*** 
(Final  Impact  and  Analysis  of  Trajectory)  and  used  in  1970-72  on  projects 
including:  1)  The  Navy  Poseidon  Reentry  Reduction,  2)  NASA  Mighty  Mouse 
Determination,  3)  Foreign  Object  Evaluation,  and  4)  System  Analysis 
Simulation. 

Early  in  1973  it  was  brought  to  SAMTEC  by  ROWS  and  implemented  on  the  360/65 
to  provide  analysis  support  for  real-time  comparisons.  Over  the  last  three 
years  it  has  been  extensively  extended  to  provide  a  reliable  yet 
versatile  reentry  program. 


*  see  reference  1  and  2 

**  see  reference  3 

***  see  reference  4  and  5 


Several  peripheral  programs  were  prepared  at  WTR  to  preprocess  information 
used  by  the  TAPP  program  as  illustrated  in  Figure  1.  The  POT  POURRI* 
program  was  written  at  WTR  to  merge  asynchronous  observation  data  from 
different  sources  onto  a  common  tape.  The  RSAT**  program  was  an  AFETR 
program  for  prefitting  atmospheric  data.  The  LRGM***  program  was  acquired 
from  Defense  Mapping  Agency  to  generate  local  a  priori  supplementary 
gravity  force  history.  The  TAL/D  program  was  an  extension  of  the  AFETR 
DMDM****  program  to  generate  high  resolution  drag  and  lift  coefficient 
history  data  from  observations.  The  POW/G  program  is  similar  to  TAL/D 
except  it  generates  thrust  and  guidance  history  for  the  boost  phase  of  the 
flight.  Finally,  the  DSTACK  and  CCPLOT  are  multiple  pass  printing  and 
plotting  programs  written  at  WTR  using  AFETR  ’destack'  philosophy. 

1.2  Characteristics 

The  TAPP  program  accepts  information  relating  the  vehicle  character, 
ambient  environmental  conditions  and  the  instrumentation  (Figure  2). 

It  provides  as  output  best  estimates  of  trajectories,  trajectory  compari¬ 
sons,  sensor  error  models,  trajectory  and  sensor  parameter  uncertainties 
and  trajectory  and  look  angle  predictions  (Figure  3). 

It  is  a  35,000  statement  FORTRAN  program  consisting  of  300  subroutines 
which  operates  in  59  overlays  on  288K  bytes  of  storage  on  an  IBM  360/65 
computer  under  IBSYS.  It  is  retained  on  tape  (and  disk)  and  manipulated 
with  a  set  of  manage  sent  routines  for  modifications  as  well  as  applications. 
Multiple  copies  of  source  libraries,  object  libraries  and  load  libraries 
are  retained  on  tape  for  backup  and  experimentation.  The  version  is 
indicated  by  the  release  date  printed  on  the  title  page  of  each  run. 

The  program  has  four  modes  of  operation  and  can  be  sequenced  (Figure  4)  by 
an  input  itinerary  of  numerical  values  2  through  5  whose  functions  are 
defined  as  follows: 


*  reference  6 

**  reference  7 

***  reference  8 

****  reference  9 
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OBSERVaTi 

Tape 


gravities 


FIGURE  2  TAPP  INPUT 


TAPP 


FIGURE  3  TAPP  OUTPUT 
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2  Parameter  Reconstruction 

3  Trajectory  Simulation/Comparison 

4  Observation  Generation 

5  Error  Analysis 

Running  times  vary  over  a  large  range  depending  on  the  problem.  A  typical 
preflight  error  analysis  with  nine  stations  and  10,000  observations 
requires  approximately  10  minutes  of  CPU  time.  A  reconstruction  run 
with  10,000  observations  may  require  from  30  to  180  CPU  minutes  depending  on 
whether  10  or  100  parameters  are  adjusted  and  upon  the  quality  of  the 
initial  conditions. 

1.3  Definitions 


The  following  is  a  partial  list  of  terms  and  their  definitions  to  provide 
the  reader  with  a  better  understanding  of  the  discussion  which  follows: 


BC 

BET 

boost 

bus 
epoch 
free  fall 

GLS 

MM 

modelled  parameter 

reconstruction 

reentry 

RV 

separation 


ball istic  camera 

best  estimate  of  trajectory 

trajectory  portion  from  lift-off  to 
final  thrust  termination 

main  body  of  launch  vehicle 
trajectory  start  point  (initial  conditions) 
portion  of  trajectory  where  no 
significant  atmosphere  is  required 

generalized  (weighted)  least  squares 

Minuteman 

parameter  which  is  subjected  to 
adjustment  in  reconstruction 

estimation,  determination,  differential 
correction  of  adjustable  parameters 

trajectory  portion  below  3Q0K  feet 
coming  down 

reentry  vehicle 

separation  of  RV  from  main  bus 
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uncertainty 


potential  parameter  error  or  measure  of 
dispersion.  Generally  the  rms  about 
the  mean. 


unmode lied 


weight 


parameter  -  parameter  not  subjected  to  adjustment 

during  reconstruction  but  whose 
uncertainty  is  propagated  into  the  final 
parameter  covariance 

certainty  (squared)  attached  to  an 
observation  or  parameter  equal  to  the 
inverse  of  the  uncertainty  squared. 
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APPLICATIONS 


The  reconstruction  logic  of  TAPP  constitutes  the  most  distinctive  capability 
of  the  program  and  has  several  range  applications. 

The  other  functions  namely  those  of  simulation  and  error  analysis  are 
primarily  subsets  of  this  reconstruction  process  with  some  extensions 
added  for  the  sake  of  convenience  and  versatility. 

2.1  Vehicle  Motion  Determination 

The  requirement  most  often  encountered  is  the  evaluation  of  a  BET,  i.e., 
best  estimate  of  trajectory,  which  is  the  motion  history  of  the  vehicle 
best  representing  available  observations.  Such  a  trajectory  is  forced  to 
be  continuous  in  position  and  velocity  by  the  using  of  equations  of  motion 
model  ling. 

For  such  a  requirement  precise  knowledge  of  the  atmospheric  profile  and 
gravitation  field  is  not  required  since  these  physical  model  variables  may 
be  treated  as  adjustable  parameters  in  the  overall  solution  (possibly  with 
some  a  priori  confidence  on  the  assumed  models  to  ensure  a  degree  of 
reasonableness) . 

For  instance,  the  atmospheric  density  versus  height  can  be  adjusted  to 
eliminate  excessive  residuals  for  any  altitude  segments  of  the  trajectory. 

An  alternative  to  this  which  would  give  the  same  overall  result  would  be 
the  estimation  of  the  drag  or  lift  coefficient  functions  versus  time  since 
they  are  both  linearly  involved  in  the  computation  of  the  fluid  forces. 

2.2  Sensor  Cal ibration 

Another  requirement  of  a  range  is  the  evaluation  of  sensor  performance 
i.e.,  checking  the  consistency  between  observations  on  a  sensor-to-sensor 
or  operation-to-operation  basis.  Simultaneous  with  the  solution  for  a 
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BET  are  the  resolution  of  individual  sensor  biases  (and  other  error  model 
terms)  to  minimize  observation  disagreement  (referred  to  as  residuals). 

Sensor  error  models  can  be  constructed  from  a  history  of  these  sensor 
error  estimates  from  several  different  BET  solutions.  Figure  5  is  a 
typical  set  of  sensor  biases  from  an  RV  separation  to  impact  solution. 

2.3  Environmental  Modelling 

The  most  prevalent  example  of  environmental  modelling  is  the  estimation  of 
gravitational  coefficients  using  multi-rev  orbital  fits.  Modelling  of 
solar  flux  and  fluid  drag  affects  for  orbital  purposes  are  also  common, 
but  the  variability  of  the  atmosphere  at  low  altitudes  severely  limits 
its  usefulness  for  reentry  or  boost  phase  environmental  modelling. 

2.4  Vehicle  Character  Estimation 


When  the  problem  is  that  of  defining  the  physical  character  of  the  vehicle 
as  well  as  its  motion,  accurate  atmospheric  observations  are  also  required. 
The  definition  of  only  the  BET  does  not  require  discerning  between  terms 
within  equation  (1)  as  to  the  specific  reason  for  the  fluid  acceleration 
magnitude  or  in  fact  even  between  the  terms  of  equation  (2). 


CD  a  p  V2 
ADRAG  '  2m 

AT0TAL  =  ADRAG  +  A,lIFT  +  ATHRUST  +  AGRAVITY 
+  ^THER 


*  where 

a 

p 

V 

m 

A 


=  drag  coefficient 

=  cross  section  area 
=  air  densi ty 

=  vehicle  velocity  relative  to  air 
=  vehicle  mass 
=  acceleration 
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When  the  physical  character  of  the  vehicle  is  to  be  determined, the  A/m 
of  (1)  is  essentially  the  required  parameter  and  knowledge  of  density  and 
fluid  velocity  (wind)  must  be  additional  observations.  Thus,  no  longer 
can  Cq  be  allowed  to  absorb  errors  in  p,  but  p  must  be  a  known  quantity. 

2.5  Predictions  -  Look  Angles,  Target  Recovery 

After  a  BET  is  formed  for  an  observed  period  it  is  sometimes  necessary  to 
extrapolate  the  trajectory  to  a  region  where  no  observations  exist  for 
purposes  of  pointing  other  instrumentation  or  even  for  target  recovery. 

The  TAPP  program  facilitates  such  runs  by  retaining  the  necessary  trajec¬ 
tory  parameters  from  a  reconstruction  run  to  apply  same  in  a  subsequent 
extended  observation  generation  run  as  depicted  in  Figure  6. 

2.6  Range  Planning 

Part  of  the  range  function  is  planning  future  efforts  including  the 
prediction  of  BET  quality  and  the  allocation  of  sensors  to  satisfy 
accuracy  requirements. 

The  basic  statistics  involved  in  a  least  squares  trajectory  reconstruction 
problem  can  be  used  to  determine  the  effects  that  specified  sources  of 
uncertainty  have  on  the  confidence  of  the  least  squares  estimated  param¬ 
eters  and  subsequently  on  the  trajectory.  These  sources  of  uncertainty, 
for  example,  may  be  systematic  inaccuracies  in  sensor  biases,  random 
variations  in  observations,  uncertainties  in  vehicle  related  parameters 
and  environmental  model  uncertainties. 

This  covariance  analysis  procedure  does  not  require  an  actual  estimation 
of  the  trajectory,  nor  does  it  require  actual  observations;  however, 
site-trajectory  geometry,  the  observation  types  and  observation  intervals 
are  required  for  each  sensor.  A  basic  covariance  analysis  of  trajectory 
uncertainties  can  then  be  obtained  by  simple  matrix  manipulation. 
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Essentially,  observation  uncertainties  are  used  to  improve  the  a  priori 
uncertainties  of  a  given  parameter  set,  which  normally  include  the  epoch 
state  vectors.  This  total  parameter  covariance  matrix  is  then  propagated 
into  trajectory  position,  velocity  and  acceleration  coordinates  for  each 
specific  time  to  provide  three  corresponding  ellipsoids  of  uncertainty  at 
each  time.  (See  Figure  7.)  Generally,  1-sigma  uncertainty  ellipsoids 
are  dealt  with  which  infer  a  20%  confidence  that  the  trajectory  is  within 
that  ellipsoid*  surface  at  that  time. 

Relative  comparisons  of  various  sensor  combinations  can  be  compared  directly 
in  terms  of  trajectory  quality  for  range  planning  purposes  although 
considerable  care  must  be  exercised  in  any  interpretation  of  such  uncertain¬ 
ties  in  an  absolute  sense  since  error  correlation  between  observations 
plays  a  very  deceptive  role  in  covariance  analyses.  For  instance,  the  use 
of  a  20  sample-per-second  (sps)  data  rate  from  a  radar  does  not  reduce  the 
vehicle  position  variance  (uncertainty  squared)  by  a  factor  of  2  from  that 
of  using  10  sps  due  to  serial  error  correlation. 


*  ”  A  2 . 8  s i gma  ellipsoid  would  have  a  95%  confidence  while  3.4  sigma 
ellipsoid  would  have  a  99%  confidence. 
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TECHNIQUES 


The  least  squares  reconstruction  process  is  basically  a  sequential  linear 
correction  process  comparable  to  the  simple  Newton-Rhapson  method  except 
in  multiple  dimensions.  Several  principal  techniques  are  discussed  belo* 
which  relate  to  its  application. 

3.1  Reconstruction  Quality  Measures 

A  strictly  quantitative  comparison  of  sensor  residuals  is  provided  by  t 
root-mean-square  magnitudes  of  each  observation  channel  as  illustrated  in 
Figure  8,  although  sensor  residual  plots  provide  the  best  qualitacive 
i nterpretation  of  reconstruction  success.  When  trended  residual  behavior 
appears  for  one  sensor  it  can  be  compared  with  that  of  another  for  trend 
correlations  which  would  likely  indicate  BET  error  rather  than  sensor 
observation  problems. 

An  example  on  Operation  7269  near  impact  {1890-1892  seconds)  showed  a  radar 
range  component  and  3  radar  angle  sets  (Figure  9)  with  highly  correlated 
non-zero  residuals.  Further,  the  range  component  of  the  ALCOR  radar 
(Figure  10)  was  also  positively  correlated  with  this  trend,  but  its 
angles  were  negatively  correlated,  i.e.,  in  contradiction.  As  it  turned 
out  multipath  had  distorted  these  ALCOR  angles  at  the  very  low  elevation 
angles  and  with  its  deweighting., the  5  sensors'  residuals  fell  neatly  to 
zero  mean  values. 

In  yet  another  case  during  reentry  a  radar  had  a  sudden  5  foot  shift  in 
range  and  disagree"  strongly  with  3  RADOTS.  In  this  case  a  sudden  vehicle 
event  caused  a  sudden  burst  of  ionization  and  shifted  the  electronic  range 
measure  leaving  the  optical  observations  unaffected. 

3.2  Observation  Weighting 

Each  observation  used  in  a  solution  can  be  assigned  a  quality  value  (an 
uncertainty,  i.e.,  actually  a  quality  inverse)  to  control  its  relative 
influence  in  the  overall  solution.  Also,  the  initial  estimate  of  each 
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FIGURE  10  RADAR  ALCOR  RESIDUALS 
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adjusted  parameter  is  assigned  such  a  quality  measure  for  the  same  reason 
(technically  each  parameter  initial  estimate  constitutes  one  added 
observation) . 

It  is  often  believed  that  any  solution  is  possible  with  the  proper 
manipulation  of  these  weights.  The  fact  is,  when  observations  are  in 
agreement,  i.e.,  consistent,  then  the  solution  is  invariant  to  weight 
changes.  That  is,  the  weighting  has  no  effect  on  the  solution.  Only  when 
observations  disagree  is  weighting  a  factor.  When  they  do  disagree  the 
cause  can  often  be  found,  using  the  residual  plots,  and  the  'erroneous1 
data  deweighted. 

In  essence  the  use  of  equal  weights,  for  disagreeing  sensors,  i.e.,  splitting 
the  difference,  is  seldon  a  satisfactory  compromise.  That  solution  is  most 
certain  to  be  a  wrong  one  since  it  agrees  with  neither  set  of  observations. 

3.3  Observation  Sufficiency 

Trajectory  solutions  constrained  by  the  equations  of  motion  require  less 
in  the  way  of  observation  constraints  at  each  time  point  than  do  point  by 
point  solutions.  For  instance,  optical  azimuth  and  elevation  observations 
over  a  sufficient  time  span  with  reasonable  geometry  can  define  a  unique 
trajectory.  If  the  geometry  js  marginal  then  a  side  condition  at  a  single 
time  point  would  normally  suffice.  As  an  example,  Figure  11  illustrates  a 
radar  tracked  bus  (with  R,  A,  E  observations)  from  which  an  optically 
tracked  RV  (A,  E  observations  only)  separates  at  some  time  prior  to  the 
observations.  Application  of  an  initial  RV  position  estimate  from  the 
bus  trajectory  constrains  the  RV  trajectory  so  that  it  then  requires  the 
adjustment  of  only  the  separation  velocity  vector  and  RV  drag  parameters. 
Varying  the  separation  time  can  further  be  used  to  minimize  RV  residuals. 

3.4  High  Resolution  Trajectory  Estimation 

The  TAPP  program  is  limited  to  estimating  25  constant  drag  parameters,  25 
horizontal  lift  parameters  and  25  normal  lift  parameters*.  Even  for 

further  limited  to  a  total  of  60  vehicle  related  parameters. 
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mildly  lifting  passive  reentry  bodies  such  as  MM  II  these  limits  yield  somewhat 
crude  results,  not  to  mention  vehicles  intentionally  maneuvered. 

A  preprocessing  program  referred  to  as  7AL/D  (Trajectory  Analysis  of  Lift 
and  Drag)  has  been  developed  which  uses  a  moving  arc  9  parameters  equations 
of  motion  filter  with  up  to  57  channels  of  observations  from  a  maximum 
of  19  sensors  to  estimate  drag  and  lift  accelerations  at  any  resolution 
permitted  by  the  data.  Although  the  position  and  velocity  vectors  are 
part  of  the  estimated  parameters  they  are  not  necessarily  useful  since  th.y 
may  include  biases  not  resolved  in  this  resolution.  A  special  precaution 
is  taken  in  this  process  to  avoid  false  acceleration  transients  caused  by 
the  starting  or  stopping  of  observation  sets  within  a  filter  span.  For 
instance,  if  radar  A  has  a  bias  and  B  does  not  then  the  apparent  trajectory 
could  appear  to  have  a  false  transient  as  depicted  by  the  dashed  line  in 
Figure  12-a.  The  constraint  is  applied  that  a  sensors  data  is  used  only 
if  it  covers  the  entire  filter  span.  With  this  constraint  the  comparable 
position  result  would  appear  as  in  Figure  12-b  and  although  position  is 
discontinuous  the  accelerations  remain  continuous.  The  situation  is 
comparable  for  rate  observations. 

The  high  resolution  values  of  C^A/W,  C^A/W  and  0  versus  time  are  recorded 
on  tape  by  TAL/D  for  suDsequent  input  into  a  TAPP  reconstruction  run 
(Figure  13-a).  Then,  in  TAPP,  a  low  frequency  adjustment  or  biasing  of 
these  profile  segments  is  accomplished  using  the  25  available  parameters 
for  both  the  drag  and  lift  components  as  depicted  in  Figure  13-b.  For 
instance,  the  net  ballistic  coefficients  for  each  time  would  be  obtained 
by  multiplying  the  values  of  profile  13-b  with  the  respective  values  of 
13-a. 
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FIGURE  13  TAL/D  A  PRIORI  BALLISTIC  COE 
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FEATURES 


Some  of  the  salient  features  of  the  TAPP  program  regarding  its  capacity, 
reliability  and  convenience  are  discussed  here.  The  more  significant 

features  are  detailed  first  and  the  remaining  are  summarized  in  the  final 
paragraph  of  this  section. 

4.1  Parameters  (modelled  and  unmodelled) 

Parameters  are  special  variables  in  the  TAPP  program  which  affect  the 
character  of  the  trajectory  and  whose  values  are  of  direct  interest  to 
the  analyst.  Such  variables  include  the  vehicle  initial  state  vectors 
and  its  drag  and  lift  character,  coefficients  relating  sensor  error  models 
and  certain  environmental  model  definitions. 

For  the  reconstruction  process  only  modelled  (adjustable)  parameters  are 
meaningful  while  in  the  error  analysis  process  unmodelled  (unadjusted) 
parameters  are  applicable.  The  use  of  unmodelled  parameters  provide  a 
means  of  including  the  uncertainties  of  reasonably  well  known  parameters, 
such  as  the  gravi tational  constant,  in  the  final  covariance  matrix  without 
subjecting  it  to  possible  alteration  in  value  caused  by  unexpected  correla¬ 
tion  problems  or  weighting  errors. 

The  available  vehicle  related  parameters  in  TAPP  are  indicated  in  Table  1 
and  the  available  sensor  and  environmental  parameters  are  shown  in  Tables  2 
and  3  respectively. 

4.2  Bounded  Least  Squares 

The  TAPP  program  does  not  rely  on  the  generalized  (weighted)  least  squares 
(GLS)  process  alone  but  utilizes  what  is  referred  to  as  a  bounded  least 
squares  process.  In  the  solution  of  a  GLS  problem  it  is  not  uncommon  for 
the  observed  weighted  residuals  to  not  follow  the  corresponding  predic¬ 
tions,  or  even  to  diverge.  Conditions  that  could  lead  to  such  circumstances 
may,  for  example,  involve  inadequacies  in  the  observational  model  or  poor 
initial  estimates  of  modelled  parameters. 


TABLE  3  SENSOR  PARAMETERS 
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Where  there  are  such  manifestations  of  nonlinearity,  it  is  necessary  to 
introduce  side  conditions  bounding  the  magnitude  of  the  correction  vector 
at  each  iteration  to  force  eventual  convergence. 

Since  the  purpose  of  the  bounds  is  to  ensure  the  good  behavior  of  each 
interation,  they  are  changed  dynamically  according  to  previous  performance. 
The  overall  weighted  sum  of  the  squares  of  all  residuals  is  this  measure 
of  iterative  performance.  If  it  grows  from  one  iteration  to  the  next  the 
solution  is  diverging,  indicating  that  the  corrections  are  too  extreme. 
Therefore,  before  proceeding,  the  last  iteration  is  repeated  with  smaller 
bounds. 

Conversely,  if  the  overall  residual  magnitude  decreases  from  one  iteration 
to  the  next  the  solution  is  converging.  Therefore,  to  allow  for  the 
possibility  of  more  rapid  convergence  through  larger  corrections  the  bounds 
are  expanded. 

To  bypass  the  bounding  process  on  any  or  all  of  the  modelled  parameters 
requires  the  use  of  negative  bound  values  for  those  respective  parameters. 

If  all  parameters  have  negative  bounds  and  an  interation  is  divergent  the 
run  will  terminate  since  no  remedial  action  can  be  taken  to  promote 
convergence. 

In  general,  the  values  of  the  bounds  should  be  set  to  that  of  the  expected 
maximum  error  in  exact  parameter.  These  bounds  are  generally  halved  when 
divergence  occurs  and  multiplied  by  1.5  when  convergence  occurs. 

Further,  discussion  of  the  bounds  logic  can  be  found  in  reference  3, 

Volume  V. 


4.3  Differential  Equations  Solution 

The  TAPP  program  solves  variational  equations  along  with  the  differential 
equations  of  motion  to  form  the  partial s  for  the  vehicle  parameters.  Two 
methods  are  available  in  TAPP  for  numerically  solving  these  differential 
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equations.  For  exoatmospheric  orbital  problems  a  10th  order  Gauss-Jackson 
technique  with  automatic  step  selection  provides  quality  solutions  in  very 
short  running  times.  For  endoatmospheric  solutions  the  differential  equations 
of  motion  are  solved  using  a  drag  linearized  Runge-Kutta  Gill  (RKG)  technique 
while  the  variational  equations  are  solved  with  only  the  RKG  method,  i.e., 
no  linearization.  Automatic  step  selection  procedures  provide  a  consistent 
solution  quality  with  respect  to  time  even  for  high  drag  parachute  drops. 

4.4  Initial  Conditions  Self  Start 

The  self-start  mode  for  deriving  initial  position  and  velocity  vectors  uses 
range,  azimuth  and  elevation  type  observations.  It  transforms  up  to  30 
triads  of  such  observations  to  Cartesian  coordinates  and  least  squares  fits 
a  0th  through  7th  order  polynomial  to  each  component.  This  method  enables 
endoatmospheric  as  well  as  exoatmospheric  reconstruction  self  starting. 

4.5  Parameter  Constraints 

Certain  adjusted  parameters  may  be  related  to  other  parameters  rather  than 
being  independent.  For  instance,  one  stations  survey  may  be  well  defined 
relative  to  another  site  which  is  itself  to  be  estimated.  In  this  case 
both  stations  are  stipulated  as  adjustable  parameters,  but  a  constraint 
matrix  is  included  to  relate  the  two  surveys.  As  a  result  the  two  surveys 
are  moved  as  a  pair  retaining  their  original  relative  locations  with  respect 
to  each  other  as  the  estimation  process  proceeds. 

4.6  Weighting  and  Editing 

All  observations  and  initial  parameter  estimates  are  weighted  to  control 
their  relative  influence  in  a  reconstruction  problem.  Observations  can  be 
weighted  as  a  group  according  to  type  and  station  or  are  individually 
weighted  on  a  point  by  point  basis  from  the  same  input  source  as  the 
observations  themselves,  i.e.,  cards  or  tape. 
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Observations  editing  can  be  controlled  based  either  upon  a  priori  input 
sigmas  or  sigmas  derived  in  the  last  iteration. 

It  should  be  noted  that  only  the  relative  magnitudes  of  the  sigmas  are 
significant  in  the  reconstruction  problem.  Their  absolute  values  are 
meaningful  only  in  the  error  analysis  problem. 

4.7  Other  Features  Summary 

Other  capabilities  of  the  TAPP  program  which  can  be  summarized  in  a  few 
words  are  listed  in  Table  4. 


TABLE  4  OTHER  FEATURE  SUMMARY 

•  MULTIPLE  VEHICLES  (SIMULTANEOUS  RECONSTRUCTION) 

•  40  SENSORS/UNLIMITED  OBSERVATIONS 

•  100  PARAMETERS  (SENSOR,  VEHICLE,  ENVIRONMENTAL) 

•  35  OBSERVATIONAL  TYPES 

•  5  ATMOSPHERE  MODELS  (DENSITY,  WIND,  PRESSURE  MODELS) 

0  1000  GEOPOTENTIAL  TERMS/ 20  POINT  MASS  MODELS 

0  MOVING  SENSOR  CAPABILITY 
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5.0 


INPUT  AND  OUTPUT 


Considerable  effort  was  invested  in  the  original  TRACE  program  and  extended 
in  the  TAPP  program  to  make  it  easy  to  use.  Large  data  sets  such  as  those 
for  the  geopotential  models  are  retained  on  tape  libraries  for  access  by  a 
single  mnemonic.  Default  constants  are  available  when  direct  input  is 
omitted.  Input  data  is  broken  down  into  various  categories  as  indicated 
in  Table  5.  The  category  input  requirement  for  each  program  function  is 
specified  in  Table  6. 

Explicit  input  error  messages  are  provided  to  clarify  mistakes  in  deck 
preparation.  Finally,  documentation  is  available  to  fully  describe  the 
methods  and  techniques  for  each  type  of  problem. 

Output  can  be  time  sequenced  or  sorted  by  station.  Special  altitude, 
latitude  or  longitude  dependent  print-outs  can  be  specified  independent 
of  time  for  trajectory  output.  Rise  and  set  information  can  also  be 
output  based  on  first  and  last  visibility  for  each  site.  Tables  7  through 
10  provide  a  summary  of  the  output  available  from  the  TAPP  program. 
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TABLE  5  TAPP  INPUT  CATEGORIES 


Name 

MODEL 

VEHICLE 

STATIONS 

SENSORS 

OBSERVATIONS 

DATA 

ATA 

COVQ 

REJECTS 

BLIST 

VI  SCON 


Constants 

Models  of  ambient  physical  conditions  which  act  as 
constraints  on  trajectory  motion  such  as  gravity 
field,  atmosphere  character,  etc. 

Vehicle  physical  descriptions  such  as  lift  and 
drag  coefficients  versus  mach,  time,  or  height, 
initial  position  and  velocity  vectors,  weights, 
areas,  thrusts,  kicks,  etc. 

Geodetic  and  astronomic  observation  or  reference 
site  positions,  ID,  indices,  etc. 

Observational  sensor  adjustable  or  fixed  parameters 
such  as  site  survey,  biases,  scale  factors,  etc. 

Observed  measurements  such  as  range,  azimuth,  ele¬ 
vation,  etc. 

Simulated  observation  generation  specifications  -  Type  II 
consists  of  the  observation  types  to  be  produced 
and  Type  I  data  generation  cards  denote  when  output 
is  to  occur. 

A  priori  (ATA)  matrix  for  modelled  parameters 
expressing  a  measure  of  their  confidence  or  quality. 

A  priori  matrix  for  unmodelled  parameters  expressing 
a  measure  of  their  quality. 

Information  specifying  observation  data  points 
or  spans  to  be  deleted  from  use  as  inputs  under 
the  OBSERVATION  identifier. 

A  matrix  of  constraints  relating  certain  modelled 
parameters  with  other  modelled  parameters. 

Visibility  contour  input  data  section. 
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TABLE  6  FUNCTION  INPUT  REQUIREMENTS 


PRINTER  PLOTS  OF  DIFFERENCES 


SIMM 


OBSERVATION  PARTIALS 


TABLE  10  OBSERVATION  GENERATION  OUTPUT 


6.0 


EXPERIENCE 


On  April  6,  1977  a  special  JR IA IG*  working  group  met  at  SAMSO  in  Los  Angeles 
to  compare  BET  results  for  the  Minuteman  III  Operation  6290  separation  to 
impact  trajectory  for  the  purpose  of  estimating  first  order  trajectory 
uncertainties**.  The  solutions  from  five  participants,  Kentron/KMR,  Lincoln 
Laboratory,  TRW  and  Xonics  were  compared  with  that  from  FEC/SAMTEC  using  the 
TAPP  program  for  a  trajectory  shape  approximately  defined  in  Figure  14.  The 
position  differences  in  terms  of  intrack  (X)  (along  earth  fixed  velocity  vector), 
crosstrack  (Y)  and  normal  (Z)  (normal  to  X  and  Y)  are  shown  In  Figure  15. 
Reference  10  details  the  parameters  of  each  participant  in  these  comparisons 
and  reference  11  details  the  6290  solution  work  accomplished  at  WTR  leading 
to  this  solution  including  comparisons  of  the  effects  using  three  different 
gravity  models  on  the  solution.  One  additional  illustration  included  at  the 
meeting  was  that  comparing  a  point' by-point  BET  with  an  equation  of  motion 
BET,  as  provided  in  Figure  16. 

In  light  of  the  success  of  this  JRIAIG  meeting,  the  1st  Strategic  Aerospace 
Division  sponsored  a  meeting  on  November  26,  1977  at  VAFB  to  compare  BET's 
for  Operation  7269,  an  MM  II  flight,  where  somewhat  fewer  observations  were 
available  from  midrange  and  uprange  sensors  (Figure  17).  Autonetics,  AVCO, 
Kentron/KMR  and  Lincoln  Lab  trajectories  were  again  compared  against  that 
generated  by  the  TAPP  program.  Position  comparisons  are  provided  in  Figure 
18.  Again  a  point  by  point  versus  equations  of  motion  comparisons  were 
illustrated  using  only  the  data  where  RADOT  and  BC  coverage  was  present 
as  in  Figure  19.  Reference  12  details  the  presentations  from  these  uncertain¬ 
ties  and  reference  13  documents  details  of  the  FEC  solution. 


*  Joint  Range  Instrumentation  Accuracy  Improvement  Group 

**  differences  from  all  possible  sources. 
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FIGURE  17  APPROXIMATE  7269  TRAJECTORY  SHAPE 
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FIGURE  18  REENTRY  TRAJECTORY  PLANE  POSITION 

DIFFERENCES  VERSUS  TIME 
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COMPARISON  WITH  KMR/KENTRON  MULTISENSOR  REPORT  POINT  BY  POINT  SOLUTION 
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BET  DEVELOPMENT  AT  WSMR 


Introduction 

The  development  of  Best  Estimate  of  Trajectory  (BET)  proqrams  at 
WSMR  began  more  than  ten  years  ago  and  these  programs  have  been  evolving 
in  several  directions.  Extensive  use  has  been  made  of  these  BET  programs 
at  WSMR  in  support  of  various  missile  and  aircraft  test  orograms.  The 
impetus  for  development  of  a  BET  program  was  provided  by  the  SAM-D 
(now  Patriot)  project  and  also  by  the  navigational  satellite  test  program. 

The  BET  program  is  still  used  in  support  of  the  Patriot  and  other  missile 
test  programs  and  has  been  used  in  support  of  several  navigational  satellite 
test  programs  including  621B,  INI  and  INHI  test  programs.  Other  test 
programs  at  WSMR  for  which  the  BET  has  been  used  to  provide  support  are 
Nike  Hercules,  the  Navy  SM2  and  ASMD  missile  test  programs,  the  Tomahawk 
and  Air  Launched  and  Cruise  Missile  (ALCM)  test  programs.  The  BET  has  also 
been  used  in  support  of  the  B1  avionics  test  program,  Lance  missile  testinc  , 
and  Completely  Integrated  Reference  Instrumentation  System  (CIRCUS)  system  tests. 

In  initial  development  of  a  BET  one  must  first  decide  as  to  what 
the  BET  must  be  able  to  do.  A  very  optimistic  but  naive  person  might 
want  his  BET  to  estimate  any  parameter  even  remotely  associated  with  a 
trajectory.  Obviously,  to  accomplish  such  a  task  would  require  a 
prohibitively  large  amount  of  instrumentation  and  computer  time.  More 
real istical ly  one  might  group  the  parameters  to  be  estimated  into  three 
classifications,  those  associated  with  the  trajectory  or  vehicle  being 
tested,  the  parameters  associated  with  the  instrumentation,  and  param¬ 
eters  associated  with  the  environment.  The  trajectory  parameters 
include  not  only  the  usual  position,  velocity,  and  acceleration  param¬ 
eters  of  the  trajectory,  but  also  parameters  such  as  attitude,  angular 
rates,  aerodynamic  coefficients,  guidance  and  control,  parameters, 
etc.  The  parameters  associated  with  the  instrumentation  include 
biases,  scale  factors,  survey  errors,  misalignments,  etc.  Atmospheric 
parameters  are  an  example  of  environmental  parameters  which  one  might 
want  to  estimate. 

From  each  of  these  groupings  one  chooses  a  relatively  small  number 
of  parameters  which  he  wants  a  BET  to  estimate.  The  actual  parameters 
which  are  estimated  on  a  given  test  might  be  a  program  option  depending 
on  the  instrumentation  available  on  the  test.  The  parameters  which  can 
be  reliably  estimated  by  a  BET  are  entirely  dependent  on  the  type  and 
geometry  of  instrumentation  available  on  the  test.  Unfortunately,  the 
high  cost  of  instrumenting  a  trajectory  and  the  relatively  few  types  of 
instrumentation  available  at  a  range  greatly  restricts  the  parameters  that 
can  be  estimated.  If  one  is  able  to  incorporate  measurements  made  on 
board  the  test  vehicle  obtained  either  by  telemetry  or  on  board  recording, 


the  possibilities  for  parameter  estimation  are  significantly  increased. 

Thus,  if  there  is  a  variation  of  instrumentation  available  from  test  to 
test  and  also  a  wide  variation  in  the  geometry  of  the  trajectory  relative 
to  the  instrumentation,  one  would  probably  want  to  include  sufficient 
program  flexibility  in  terms  of  options  on  the  parameters  to  be  estimated 
to  match  the  variable  instrumentation.  In  the  early  development  of  a  BET 
at  WSMR  there  was  considerable  interest  in  developing  a  highly  flexible 
program  that  could  utilize  Patriot  telemetry  in  conjunction  with  range 
instrumentation  data  to  do  a  relatively  complete  estimation  of  aerodynamic 
and  control  parameters  of  this  system.  To  develop  a  BET  in  this  way 
would  require  an  extreme  amount  of  program  flexibility  if  the  program  were 
to  be  used  for  other  than  Patriot  trajectories.  Thus,  we  decided  to  pur¬ 
sue  a  more  general  development  of  a  BET  that  required  much  less  flexibility 
and  would  apply  to  the  wide  variety  of  trajectories  flown  at  WSMR  and  be 
able  to  utilize  any  available  range  instrumentation  and  have  a  limited 
amount  of  ability  to  utilize  measurements  made  on  board  the  test  vehicle. 

We  currently  have  two  operational  BET  programs  at  WSMR.  The  first 
of  these  two  programs  to  be  developed  was  our  Kalman  filtering  3ET.  This 
program  has  been  operational  since  1971  and  has  been  used  extensively  in 
support  of  WSMR  test  programs.  It  is  seldomly  used  at  the  present 
time  because  we  currently  favor  the  use  of  our  other  BET  program  for 
reasons  which  we  will  discuss.  Our  second  BET  program  is  a  nonlinear 
least  square  batch  processing  type  of  algorithm  which  we  shall  call 
MISTE.  We  began  development  on  the  MISTE  program  about  1973  and  it 
became  operational  in  1975.  Since  m id— 1 97 6  it  has  been  used  almost 
exclusively  for  BET  support  at  WSMR. 

The  intention  of  this  paper  is  to  provide  an  overview  description  of 
the  WSMR  BET  programs  including  coordinate  systems  used,  trajectory 
models,  basic  flow  diagrams,  basic  estimation  equations,  input  require¬ 
ments,  output  options,  data  editing,  etc.  Our  BET  methods  are  under 
continuous  development.  One  of  the  most  important  aspects  of  our  BET 
efforts  is  continuous  improvement.  Thus,  this  paper  will  also  describe 
our  current  directions  of  efforts  at  improving  our  BET  programs.  Another 
important  aspect  of  WSMR  BET  efforts  has  been  programs  which  were  developed 
as  derivatives  of  the  BET  programs.  These  developments  will  also  be  described. 

Coordinate  Systems 

Many  different  coordinate  systems  are  used  in  BET  data  reduction. 

The  two  most  prominent  coordinate  reference  systems  in  the  WSMR  BET's 
are  the  WSCS  and  the  BET  systems. 
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The  WSCS  coordinate  system  has  its  origin  at  33°  05'  latitude, 

106°  20'  W  longitude  and  altitude  at  sea  level.  The  x-axis  is  positive 
east,  y-axis  positive  north,  and  the  z-axis  positive  up.  The  WSCS  system 
is  the  basic  cartesian  system  used  at  WSMR  and  all  instrument  sites, 
launch  sites,  etc.,  are  referenced  to  WSCS  coordinates. 

The  trajectory  coordinates  estimated  by  the  BET  program  are  referenced 
to  the  BET  cartesian  system.  The  BET  system  is  the  basic  coordinate 
system  of  the  BET  program.  The  origin  of  the  BET  system  is  at  latitude 
Ag,  longitude  yg,  and  altitude  Hg.  The  x-axis  is  positive  east,  y-ax.s 

positive  north,  and  z-axis  positive  up.  Equivalently  the  origin  of  the 
BET  system  may  be  defined  by  specifying  its  WSCS  coordinates  (Xg/W, 

yB/w»  z b/W ) ‘  For  a  9round  launched  missile  trajectory  the  most 

commonly  chosen  origin  for  the  BET  is  the  launch  point.  For  an  aircraft 
trajectory  it  is  common  to  choose  the  BET  system  to  be  the  WSCS  system. 

The  matrix  which  rotates  from  the  BET  system  to  the  WSCS  system  will  be 
denoted  by  MWg.  Mwg  is  computed  from  a  knowledge  of  the  latitudes  and 

longitudes  of  the  BET  and  WSCS  systems.  The  specific  form  of  this  3x3 
rotation  matrix  will  not  be  given  here  but  may  be  found  in  [21]. 

The  position  translation  vector  from  the  BET  system  to  the  WSCS  system 
has  components  (xg^,  yB/w>  zb/^)  and  are  the  components  of  the  BET 

origin  in  the  WSCS  system. 

Let  (x,  y,  z)  be  the  velocity  components  in  the  BET  system  of  the 
target  reference  point  with  respect  to  the  BET  origin.  Let  the  origin 
of  a  trajectory  Cartesian  system  be  located  at  some  reference  point  on  the 
target.  Define  the  rotation  matrix  from  the  BET  system  to  the 
trajectory  system  as 


-  - 

—  — 

V 

x 

0 

=  mtb 

y 

0 

• 

z 

- 

_ 

where  v  = 


(x2  +  y2  + 


z2) 


i 
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The  unit  base  vector  tj  of  this  trajectory  system  is  tangent  to  the 

target  trajectory.  The  orientation  of  the  trajectory  system  base 
vectors  along  a  trajectory  is  shown  below. 


Note  that  the  transformation 


is  singular  when  Vg= 0, 


i . e. ,  when  the 


trajectory  is  vertical.  In  order  to  avoid  this  singularity  we  choose  a 
different  trajectory  coordinate  system  for  nearly  vertical  trajectories. 


At  each  range  instrumentation  site  we  define  a  cartesian  coordinate 
system  with  origin  having  the  geodetic  coordinates  (x.,  v  . ,  H.)  and  with 
the  x-axis  positive  east,  the  y-axis  positive  north, 
and  the  z-axis  positive  up.  The  rotation  matrix  from  an  instrument 
coordinate  system  to  the  BET  system  is  denoted  by  MD.  and  is  computed 

d! 


from  the  latitude  and  longitude  of  the  two  origins.  The  coordinates 
of  the  origin  of  an  instrument  coordinate  system  with  respect  to  the  BET 
system  origin,  denoted  by  (x^g,  y.^g,  z.^g),  are  computed  from 
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Xi/B 

Xi/W 

XB/W 

yi/B 

=  M  ' 

nBW 

yi/W 

- 

yB/W 

Zi/B 

Zi/W 

ZB/W 

where  (x^,  y^w>  z.^)  are  the  coordinates  of  the  origin  of  the 
instrument  coordinate  system  with  respect  to  the  WSCS  origin. 


The  trajectory  estimated  by  the  BET  program  must  refer  to  some 
specific  point  on  the  target.  In  a  single  instrumentation  system  type 
of  data  reduction  this  target  reference  point  would  usually  be  the  point 
on  the  target  to  which  the  measurement  is  referenced,  e.g.,  DOVAP  antenna, 
radar  antenna,  painted  cross  for  optics.  For  a  BET  with  several  types  of 
measurement  systems  tracking  the  target  it  is  absolutely  necessary  to 
translate  each  measurement  to  a  common  point  on  the  target.  This  point 
may  be  specified  by  the  range  user.  This  point  is  called  the  target 
reference  point  and  is  the  origin  of  the  target  coordinate  system.  For 
an  elongated  target  such  as  a  missile  or  an  aircraft  the  target  1-axis  or 
s-axis  is  parallel  to  the  long  axis  and  positive  toward  the  nose.  The 
target  2-axis  or  w-axis  is  normal  to  the  s-axis  and  positive  upward.  The 
3-axis  or  b-axis  completes  a  right  handed  system. 


Since  the  target  coordinate  system  is  rigidly  tied  to  the  target,  the 
orientation  of  this  coordinate  system  is  unknown  unless  there  are  attitude 
or  gimbal  angle  measurements  relating  the  target  orientation  to  a  coor¬ 
dinate  system  of  known  orientation.  In  the  absence  of  such  measurements, 
it  is  necessary  to  assume  that  the  target  coordinate  system  axes  are 
parallel  to  the  trajectory  coordinate  axes  previously  defined.  If  gimbal 
angle  measurements  are  available,  the  transformation  from  the  target 
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coordinate  system  to  the  BET  system  can  be  computed.  The  computation  of 
this  transformation  is  peculiar  to  the  type  of  inertial  measurement  unit 
( IMU )  being  used.  If  no  angle  measurements  are  available,  it  is  necessary 

to  set  =  Nljg.  If  the  coordinates  of  a  measurement  reference  point  in 

the  target  coordinate  system  are  (d_,  d  ,  d,  ),  then  the  coordinates  of  thi 

S  W  D 

point  in  the  BET  system  are 


1 

—  ") 

d 

d 

x 

s 

d 

y 

MBA 

d 

w  1 

d 

d.  ; 

z 

b 

Measurement  Models 


Let  (x,  y,  z)  be  the  coordinates  of  the  target  reference  point  in  the 

BET  coordinate  system.  Also,  let  (d  d  .,  d  .)  be  the  BET  coordinates 

xi  yi  zi 

of  the  measurement  reference  point  on  the  target  (radar  antenna,  optical 
cross,  etc.)  with  respect  to  the  origin  of  the  target  coordinate  system. 


d  ■  1 

xi 


d  . 
yi 


Z1 


=  M 


BA 


r 

!  d  . 

I  ST 


d  . 
wi 


dbi 


(5) 


where  is  the  rotation  matrix  from  the  target  coordinate  system  to 

the  BET  coordinate  system  and  (d  . ,  d  d.  .)  are  the  coordinates  of  the 

measurement  reference  point  in  the  target  coordinate  system.  Let  the 
coordinates  of  the  measurement  reference  point  on  the  target  with  respect 
to  the  instrument  coordinate  system  be  denoted  by  (x-jy . ,  y^  . ,  z^.). 
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They  are  given  by 


- 

XT/i 

*  +  dx1  -  Xi/B 

yT/i 

=  MB1 

y  *  dyi  '  *1/B 

ZT/i 

z  +  d2i  '  Z1/B 

L 

where  M^.  is  the  rotation  matrix  from  the  local  instrument  coordinate 
system  to  the  BET  coordinate  system  and  (x./D,  y./D,  z./D)  are  the 

l/D  1/d  1/d 

coordinates  in  the  BET  system  of  instrument  origin  given  by  (3). 

Radar  Measurements 

The  range,  azimuth,  and  elevation  measurements  of  the  FPS-16  and 
MPS-36  radars  are  ideally  modelled  as 


R  -  (xT/i 

2  2  ^ 

+yT/1  *  zT/i> 

(7) 

A  =  tan1 

XT/i 

yT/i 

(8) 

E  =  tan1 

ZT/i 

(9) 

— 2 j- 

XT/i  +  yT/i^ 

These  are  ideal  models  of  the  measurements  which  assume  that  the 
systematic  errors  such  as  refraction,  zero  set,  collimation,  tilt, 
etc.  have  been  removed  from  the  measurements  by  calibration.  Although 
the  radar  measurements  have  been  calibrated  and  the  estimated  systematic 
errors  removed  from  the  data,  there  is  usually  a  significant  amount  of 
systematic  error  remaining  in  the  data.  Thus,  one  must  use  an  error 
model  of  the  radar  measurement  and  estimate  the  unknown  parameters  of 
this  error  model.  A  fairly  complete  error  model  for  radar  measurements  is 
given  by 


AR  =  r0  +  rift 


(10) 
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aA  =  a0  -  a3tanEcosA  -  a2tanEsinA  +  a3secE  +  a4A  (11) 

aE  =  e0  +  a^inA  -  a2cosA  +  e3E  (12) 


In  the  above  error  model  (r0,  a0,  e0)  are  zero  set  error  coefficients, 

a3  and  a2  are  tilt  errors,  a3  is  a  collimation  error,  and  r3,  a4,  ei 

are  a  combination  of  timing  error  and  servo  lag  errors.  Other  error 
terms  could  also  be  included  in  this  model.  The  relative  geometry  of 
the  trajectory  and  instrumentation  plan  will  determine  which  coefficients 
in  the  error  model  can  be  reliably  estimated.  For  most  trajectories  at 
WSMR,  it  is  not  possible  to  estimate  more  than  just  the  zero  set  error 
coefficients  (ra,  aG,  e0).  This  situation  is  acceptable  if  the  other 
error  coefficients  have  L  ■‘n  reliably  estimated  by  the  prsmission  cali¬ 
bration  procedure.  For  many  WSMR  trajectories  the  relative  geometry  will 
not  even  allow  the  estimation  of  the  zero  set  coefficients.  Some  of 
the  radars  available  at  WSMR  have  the  capability  to  provide  a  range  rate 
measurement.  The  ideal  model  for  this  measurement  is 


ft= 


xT/ixT/i 


yT/iyT/i 


zT ,  .zT  . . 
.T/  0/1 


(13) 


No  errorsare  estimated  for  range  rate  measurements.  However,  an 
approximation  for  range  rate  refraction  is  removed  from  the  data. 


Optical  Measurements 


The  fixed  cameras  and  cinetheodol ites  provide  measurements  of  the 
azimuth  and  elevation  angles  of  the  line  of  sight  from  the  camera 
coordinate  system  origin  to  a  point  on  the  target.  The  ideal  measurement 
models  are 


A  =  tan1 

yT/i 


E  =  tan1 


-T/i. 

,  2  2  > 

xT/i  +  yT/i) 


i 


(H) 


(15) 
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The  above  model  assumes  that  systematic  components  of  error  have  been 
reliably  estimated  by  premission  calibrations  and  removed  from  the  data 
and  that  the  measurements  have  been  corrected  for  refraction.  The  pre¬ 
mission  optical  calibrations  are  usually  reliable  so  that  it  is  not 
necessary  to  estimate  coefficients  in  a  complete  error  model.  At  most 
it  may  be  necessary  to  estimate  zero  set  error  coefficients  from  the 
trajectory  measurement  data. 


Measurements 


The  Dovap  measurement  system  is  a  two  way  CW  doppler  system  which 
measures  a  loop  range  change  or  alternatively  an  average  loop  range 
rate  of  a  target  between  consecutive  sampling  times  t^  and  t^+1> 

Although  the  Dovap  system  has  been  phased  out  at  WSMR,  it  has  been  used 
extensively  in  BET  reductions  at  WSMR  and  well  illustrates  the  use  of 
doppler  measurements  in  a  BET. 


Let  (Xy,  y-j.,  Zy)  be  the  coordinates  of  the  ground  transmitter 
antenna  in  the  BET  coordinate  system.  Also,  let  (xR,  yR,  zR)  be  the 
coordinates  of  a  ground  receiver  antenna.  Let  (d^y,  d^y,  d^y)  and 
(dxR,  dyR,  dzR)  be  the  BET  coordinates  of  the  Dovap  transmitting  and 
receiving  antennas  on  the  target.  Let  (dsy,  dwy,  dfay)  and  (d$R,  dwR,  dbR) 
be  the  coordinates  of  these  antennas  in  the  target  coordinate  system. 

Then 


’13 


(18) 


The  loop  range  change  model  for  a  Dovap  receiver  measurement  is 
s(tr  tj+1>  .  RT(t1+1)  ♦  RR(tj+1)  -  RT(tj)  -  RR(tj) 


where  Rj(t^)  is  the  range  between  the  ground  and  target  transmitter  antennas 
at  time  t^  and  RR(t..)  is  the  range  between  the  ground  and  target  receiver 
antennas  at  time  t^ . 


+<z(ti>*w-zR)2:ii 


(19) 


RT ( ti )  *  [(x(t1)+dxT(tf)-xT)2t(y(t1)+dyT(tj)-yT)2 
+(x(tj)+dzT(tj)-2T)2]J 


(20) 


The  Dovap  measurements  are  corrected  for  refraction  as  described  in  [21]. 
No  other  systematic  errors  are  assumed  to  be  present  in  the  Dovap 
measurements. 

An  alternative  to  the  above  position  type  measurement  model  for  Dovap 
is  a  velocity  type  measurement  model.  If  the  dovap  measurements  are 
divided  by  (t^+j-  t . ) ,  the  result  is  an  average  loop  range  rate  over  the 

interval  (t. ,  t^+]).  This  average  loop  range  rate  can  be  used  to  approxi¬ 
mate  the  loop  range  rate  at  tAVE=  (t..+  t.+i)/2.  The  loop  range  rate  model 
of  the  Dovap  measurements  is 

£  =  (x^xT)(x-HdxT-xT)+(y+jyT)(y+dyT-yT)  +  (i+3zT)(z^dzT-zT) 

RT  (21) 

+  (x+dxR)(x^dxR-xR)  +  (^vR)(y+dyR-yR)  +  (z^zR)(z^dzR-zR) 

RR 


214 


In  the  above  model  all  time  varying  quantities  are  evaluated  at  t^£. 

The  quantities  (d'x^,  d^)  and  ( d^y ,  d*^,  d"^)  can  be  approximated 

by  differencing  or  differencing  and  smoothing  the  quantities  defined 
in  (16)  and  (17). 

LTN-51  Inertial  Measurement  Unit 

The  LTN-51  inertial  measurement  unit  (IMU)  is  one  example  of  an  1MU 
whose  measurements  have  been  used  in  a  WSMR  BET.  The  LTN-51  provides 
inertial  acceleration  measurements  from  accelerometers  mounted  on  a 
4-gimbal  stable  platform.  In  addition,  the  angle  readings,  which 
measure  the  orientation  of  the  target  with  respect  to  a  coordinate  system 
established  by  the  platform,  are  also  used  in  the  BET  program.  In  terms 
of  the  tangential,  normal  and  lateral  accelerations  along  the  trajectory 
the  measured  accelerations  are  (approximately) 


where  (A-p  A^,  A^)  are  the  acceleration  components  in  the  trajectory 
coordinate  system,  (d, ,  d„,  d3)  are  constant  unknown  scale  factor  errors 
which  are  to  be  estimated,  (&l,  62,  6^)  are  constant,  unknown  platform 
misalignment  errors  to  be  estimated,  and  (bp  b2,  b3)  are  constant  bias 
errors  to  be  estimated.  The  3x3  matrices  MD  T,  D,,  and  S, ,  are  known 

r  q  I  U  U 

quantities  depending  mainly  on  the  latitude  and  longitude  of  the  platform, 
the  latitude  and  longitude  of  the  origin  of  the  BET  coordinate  system, 
and  the  velocity  of  the  target. 


Other  Measurements 

Other  types  of  measurements  such  as  telemetered  acceleration,  angle 
measuring  equipment  (AME)  measurements  which  are  direction  cosines,  velocim- 
eter,  and  laser  tracker  (R,  A,  E)  measurements  have  been  used  at  WSMR  in  the 
BET.  The  laser  tracking  measurements  will  probably  be  more  often  used  in  the 
future.  It  is  a  relatively  simple  matter  to  add  capability  to  either  of  the 
BET  programs  so  that  they  can  process  additional  types  of  measurements. 
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Kalman  Filter  BET 


The  Kalman  filter  assumes  a  dynamic  model  of  the  trajectory  in  terms 
of  a  state  vector  x.  A  model  of  the  systematic  measurement  error  compo¬ 
nents  is  specified  as  a  function  of  a  bias  vector  b.  Also,  a  measurement  model 
is  specified  as  function  of  the  state  vector  x  and  the  bias  vector  b.  The 
state  vector  x  is  takento  be  a  9-vector  of  positions,  velocities,  and 
accelerations.  The  state  vector  and  trajectory  model  are  given  by 
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(23) 


2  2  1  o  1 

where  v^  =  (x^  +  x ).')“,  v  =  (vg  +  x^)J  and  Ay,  A^,  AL  are  the  components  of 

acceleration  in  the  trajectory  coordinate  system.  Thus,  we  are  modelling 
the  trajectory  as  a  constant  acceleration  system.  We  compensate  for  this 
mismodelling  by  adding  an  uncertainty  term  to  the  state  equation  so  that 
we  actually  assume 


x  =  f(x)  +  w 


(24) 


where  w  is  a  random  vector  with  zero  mean  and  covariance  matrix  Q. 

Q  is  a  9x9  matrix  with  zero  entries  except  for  the  7th,  8th,  and  9th 
diagonal  entries  which  we  use  to  compensate  for  mismodelling  the  accel¬ 
eration. 
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The  Kalman  filter  provides  an  optimal  estimate  of  the  state 
x(tk)  of  a  linear  system  corrupted  by  additive  Gaussian  noise  con¬ 
ditioned  on  all  observations  up  to  time  t^.  In  trajectory  estimation 

we  do  not  have  a  linear  system  but  rather  a  nonlinear  system  due  to 
the  nonlinear  nature  of  the  trajectory  equations  and  the  nonlinear 
relationship  between  the  measurements  and  the  state  variables. 

The  Kalman  filter  equations  can  be  extended  to  this  nonlinear 
situation  and  while  the  resulting  equations  and  estimates  are  no 
longer  optimal,  they  still  provide  good  estimates  which  are  easily 
implemented  on  the  computer.  The  resulting  filter,  which  is  called 
an  extended  Kalman  filter,  is  a  most  popular  way  of  implementing  a 
recursive  estimation  procedure  for  nonlinear  systems  and  has  been  the 
subject  for  considerable  development  and  evaluation  in  engineering 
literature,  [1]  -  [4], 

Let  x(tk)  be  the  state  vector  of  the  trajectory  at  time  t^.  The 
trajectory  dynamics  are  governed  by  the  nonlinear  dynamic  equation 

x(tk)  =  f<x(tkj)  +  wk  (25) 

where  f(«)  is  a  known  nonlinear  function  specified  in  (23)  and  wk  is  a 
vector  of  Gaussian  noise  with  mean  zero  and  covariance  matrix  Qk<  The 
ith  trajectory  observation  at  time  tk  is 

ZA}  =  hi(x(V) +  gi(tk)b(tk}  +  vi(tk5  (26) 

where  h.(*)  is  a  known  nonlinear  function,  v.(t. )  is  zero  mean  Gaussian 

1  I  K 

noise  with  covariance  R. (t^),  b( t^ )  is  a  constant  m-vector  of  biases  to 

be  estimated,  and  gl(tk)  is  a  known  lxp  matrix.  We  assume  and  process 

only  scalar  observations  with  no  loss  of  generality.  We  use  the  decom¬ 
position  procedure  of  B.  Friedland  [5]  -  [7]  to  separate  state  estimation 
from  bias  estimation. 


In  the  following  we  denote  by  x(k)  the  estimate  at  time  t,  given 

A  K 

all  measurements  up  to  time  t.  and  by  x(k/k-i)  the  estimate  at  time  t 

K  -  k 

given  all  measurements  up  to  time  t^.  b(k)  is  the  estimate  of  the 

bias  vector  given  measurements  up  to  time  tk . 
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Thus,  x(kjk-l)  is  an  estimate  predicted  from  x(k-l).  x(k)  is  an  estimate 
obtained  assuming  no  biases  are  present.  The  estimation  equations  for 
our  Kalman  filter  BET  are 

At  a  Time  Update 


x(k+l  |k)  =  x(k)  +  f(x(k))Atk+]  +  J(x(k))f(7(k))— jp1 

J(^))=  C^x(k) 


(27) 


Pk+l/k  \^Pk  +  •5QkAtk+lK+  ,5Qk;tk+l 
\  =  I  +  J(x(k))Atk+1  +  J2(x(k))  — p 


(28) 


(29) 


In  the  above  Px  is  the  9x9  covariance  matrix  of  x(k)  and  P^k  ^  is  tne 
covariance  matrix  of  x(k/k-l).  In  obtaining  the  prediction  equation 


(27)  for  x  a  second  order  Taylor  series  integration  of  (23)  was  used. 

The  matrix  Ricatti  equation  which  governs  the  evolution  of  the  covariance 
matrix  is  integrated  using  a  trapezoidal  rule  to  obtain  (28). 


Zero  Bias  Filter  at  Measurement  Update 
J1)  Ji-lp  T 
Pk  '  ‘he  * 


A 


i  =  l  ,m 


(30) 


JO  JOt  (i-Or 

w(k)  =  Pk  Hj/(Ri(tk)  +  HiPk  H]) 


>  m=#  measurements  at  t, 

I  k(31) 


Jl)  Ji-1)  JO  Ji-1) 

x(k)  =  x(k)  +  w(k)  (zi(tk)  -  h.(x(tk))) 


Hi 


Ji-0 

x=x(tk) 


JO)  _ 

X  =  x(k/k-l) 


J 


(32) 
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Bias  Filter  at  Measurement  Update 


-  !>X>  +  HlWk> 


a1-(k)  =  R.j(tk)  +  HiPk^1  iJhT  =  variance  of  residual 


P£i}(k)  =  (Pji_l)(k)  +  sis  ./a • (k) ) 


i  r  i 


h= i 


(i) 


(k)  =  pji_l)(k)s1(k)/(a1(k)  +  sjp(1->(k,s1, 


_( i - 1 ) 

^(k)  -  z i ( k )  -  h(x(k)) 


=  residual  from  zero  bias  filter 


b(l)(k)  =  b(w)(k)  +  w^l}(k)  (r.(k)  -  sT(k)b(i_l)(k)) 


>m 


COMBINING  MATRIX  UPDATE 


:o) 


T i ( k )  =  Ti_1(k)  -  w  (k)  sl(k)  i=i,m 


(33) 

(34) 

(35) 

(36) 

(37) 

(38) 


(39) 


OPTIMAL  ESTIMATE 

x(l)(k)  =  x(l\k)  +  T.(k)b(l)(k) 


i  =  i  ,m 


(40) 


The  direct  implementation  of  the  above  Kalman  filter  equations  will 
often  lead  to  serious  numerical  stability  problems.  This  instability, 
called  filter  divergence  in  the  literature  [8]  -  [10],  is  primarily 
caused  by  mismodelling  of  the  target  dynamics,  mismodelling  of  the 
measurements,  or  by  errors  in  the  noise  statistics  used  by  the  filter. 

The  divergence  due  to  these  sources  have  been  successfully  treated  by 
several  methods,  most  commonly  by  the  addition  of  ficticious  process 
noise  covariance  Another  contributor  to  the  filter  divergence  problem 
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is  numerical  inaccuracies  in  the  filter  computations.  These 
numerical  inaccuracies  may  act  in  such  a  way  that  the  filter  covariance 
matrices  lose  their  positive  definite  character  causing  divergence  of 
the  filter.  One  very  effective  method  of  preventing  divergence  caused 
by  numerical  errors  is  to  implement  a  square  root  version  of  the 
Kalman  filter  equations  [5],  [6],  [11]  -  [14].  A  square  root  filter 
computes  a  lower  triangular  square  root  matrix  L  such  that 

P=LLT  or  such  that  P=LDL^  where  D  is  a  diagonal  matrix  and  Pisa 
covariance  matrix.  In  the  square  root  implementation  the  square  root 
matrix  is  computed  and  propagated  in  time  rather  than  the  covariance 
matrix.  All  filter  computations  are  expressed  in  terms  of  the  square 
root  matrix.  Our  Kalman  filter  BET  is  implemented  by  using  a  numerically 
efficient  square  root  formulation  of  the  Kalman  filter  equations.  The 
result  is  a  stable  filter  that  is  almost  as  numerically  efficient  as 
the  basic  Kalman  filter  computations  but  has  much  more  accuracy. 

The  dimension  of  the  bias  vector  b  can  be  quite  large  sometimes  as 
large  as  40.  Then  the  computation  of  and  computations  involving  the 
bias  covariance  matrix,  for  each  scalar  measurement  and  each  time  point 
along  the  trajectory  sometimes  presents  a  great  computational  burden. 

The  implementation  of  the  Kalman  filter  for  a  BET  is  certainly  not 
difficult.  There  are,  however,  two  important  features  which  must  be 
implemented  in  any  BET  program  and  whose  successful  implementation  is 
a  must  for  good  BET  performance.  These  two  features  are  measurement 
data  editing  and  adaptive  filtering.  In  our  Kalman  filter  BET  and  in 
other  Kalman  filtering  programs  a  rather  simple  method  of  measurement 
editing  has  been  used.  In  this  method  an  observation  whose  absolute 
predicted  residual. 


_( i - 1 )  T 

jzi(t|<)  -  H.x(k)  -  g ■  (k)b(k) ] 
is  greater  than, 

(i-Dr 

aa^k)  =  a(Ri(t|<)  -  HiPk  H* ) , 

where  a  is  a  suitable  scalar  usually  in  the  range  3<a<6,  is  not  processed. 
This  is  an  oversimplification  of  the  measurement  editing  process  which 
at  least  in  a  multi-instrument  BET  Kalman  filter  must  be  considerably 
more  complex.  In  a  Kalman  filtering  BET  the  editing  must  attempt  to 
answer  the  questions:  Did  the  residual  fail  the  criterion  because  of  a 
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jump  in  bias,  because  of  a  sudden  change  in  noise  level,  or  because 
of  a  localized  wild  data  point?  How  many  consecutive  times  may  a 
measurement  fail  the  criterion  before  we  decide  the  failures  are  not 
localized  wild  data  but  are  caused  by  a  bias  jump  or  noise  level  jump? 
What  do  we  do  if  we  decide  a  bias  or  noise  jump?  Should  measurement 
noise  covariance  be  modified  when  criterion  failure  occurs?  Our 
measurement  editing  technique  answers  all  these  questions  and  while 
quite  successful,  it  consists  of  several  ad  hoc  techniques  and  special 
conditions  peculiar  to  measurement  type.  A  more  unified  method  of 
measurement  editing  would  be  desirable.  We  are  now  working  on  some 
filtering  techniques  which  will  be  discussed  later  where  measurement 
editing  is  inherent  in  the  filter  derivation.  Except  for  some  special 
conditions  steps  in  the  BET  measurement  editing  are  given  by  the 
fol lowing: 


_  T  a/4  l) 

(1)  If  | ri ( k)  -  g.j Ck ) b v  '(k)|  £4or,  orocess  measurement  and  update 


computation  of  R^t^), 

(2)  If  4or<Iri  (k)  -  g{(k)bVl'°(k)|112jr, 
update  computation  of  R^(tk), 


do  not  process  measurement. 


but 


(3)  If  |r^(k)  -  gj(k)b^-1)(k)  |>12°r,  do  not  process  measurement  and  do 
not  update  computation  of  R-j  ( t ^ ) » 


In  the  above,  is  the  variance  of  the  residual,  r-j(k) 


T  (i-i) 

g i ( k ) b  (k). 


(4)  If  condition  in  (2)  occurs  eight  consecutive  times,  decide  bias  has 
changed  and  reinitialize  corresponding  zero  set  bias  compo  nt  in  bias 
filter  with  old  bias  estimate  plus  average  of  last  eight  residuals. 

(5)  If  condition  (3)  occurs  four  consecutive  times,  assume  instrument 
failure  and  reinitialize  filter  to  reflect  this  condition.  Measurement 
may  try  to  reenter  solution. 

The  second  filter  feature  which  is  necessary  for  good  performance 
is  that  the  filter  must  be  able  to  change  its  parameters  in  order  to 
adapt  to  changes  in  measurement  noise,  process  noise,  and  errors  made 
in  mismodelling  dynamics  and  measurements.  Compensation  for  errors  in 
modelling  are  usually  made  thru  proper  selection  of  the  process 
covariance  Q^.  Since  the  principle  source  of  modelling  error  in  the 

WSMR  BET  filter  is  due  to  the  assumption  of  constant  accelerations ,  we 
add  process  noise  only  to  these  acceleration  states  by  choosing  the 
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7th,  8th,  and  9th  diagonal  elements  of  the  matrix  Q^.  All  other 

elements  of  Qk  are  zero.  We  choose  the  level  of  these  elements 

depending  upon  the  level  of  acceleration  (assuming  high  levels  must 
mean  a  high  future  rate  of  change  of  acceleration)  along  various  parts 
of  the  trajectory.  This  is  not  done  automatically  in  our  filter,  but 
is  set  either  by  a  knowing  what  acceleration  levels  are  expected  for 
the  mission  or  by  running  a  preliminary  filter  over  the  trajectory 
data  to  determine  approximate  acceleration  levels.  The  measurement 
noise  covariance  used  by  the  Kalman  filter  program  is  computed  for 
each  measurement  by  applying  a  fading  memory  filter  to  the  scalar 
measurement  residuals.  There  have  been  many  efforts  to  construct 
Kalman  filter  algorithms  which  will  automatically  adapt  to  maneuvering 
targets.  However,  none  of  the  methods  we  have  seen  seem  to  be  satis¬ 
factory  for  our  Kalman  filter.  There  are  some  current  research  efforts 
in  adaptive  estimation  which  offer  sens  promise  and  we  are  eagerly 
awaiting  the  results  of  the  research. 

Batch  Proic  ;or  BET 

Let  x ( t, )  be  a  state  vector  for  a  trajectory.  If  we  have  only 
position  measurements  available,  we  let  x(tk)  be  the  3-vector  with 

components  (x,  y,  z),  i.e.,  the  BET  coordinates  of  the  trajectory. 

If  in  addition  to  position  measurements,  we  also  have  range  rate 
measurements  available,  we  let  x(t^)  be  a  6-vector  with  components 

(x,  y,  z,  X,  y,  2).  The  batch  processor,  which  we  call  MISTE, 
combines  measurements  from  radars ,  fixed  cameras,  cinetheodol i tes ,  and 
laser  trackers  to  estimate  the  state  vector  x(tk)  of  the  target. 

Attitude  measurements  may  also  be  used  in  the  batch  processor  to  relate 
each  measurement  to  the  common  target  reference  point.  Batch  processor 
in  this  case  means  that  measurements  for  the  entire  trajectory  are 
processed  simultaneously  rather  than  sequentially  in  time.  We  obtain 
an  estimate  x(tk)  simultaneously  for  each  t^  along  the  trajectory  and 

also  obtain  estimates  of  the  measurement  bias  terms.  Acceleration 
states  or  velocity  and  acceleration  states  of  the  trajectory  are  obtained 
by  running  an  adaptive  variable  lag  smoother  over  the  raw  trajectory 
estimates  obtained  from  the  batch  processor. 

There  are  some  advantages  to  the  batch  processor  formulation  of  the 
BET.  Since  no  trajectory  model  is  assumed  by  the  batch  processor,  there 
are  no  errors  in  the  position  estimates  or  in  the  measurement  bias 
estimates  induced  by  trajectory  mismodel 1 ing.  The  influence  of  mismodelling 
errors  on  position  and  bias  estimates  are  very  significant  in  the  Kalman 
filter  BET.  We  also  have  found  that  the  batch  processor  BET  is  considerably 
more  efficient  in  terms  of  computer  time  than  the  Kalman  filter  BET. 
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(41) 


/ 


Let  the  scalar  measurement  model  be  as  before,  i.e., 
zi Ctk)  =  Cx  Ctk) )  +  g^(k)b  +  v(tk) 

Assume  we  have  measurements  along  the  trajectory  at  times  tk>  k=l ,  2, —  ,n 
and  for  i=l,2, —  ,m.  Again  let  R.j(tk)  be  the  variance  of  v(tk).  Then  the 
batch  processor  estimates  the  p-vector  b  and  the  state  vector  x(tk)  so  that 


(z.(tk)  -  h-(x(tk)  -  gl(k)b) 


(42) 


is  minimized.  Since  we  have  a  nonlinear  estimation  problem  due  to  the 
nonlinear  relation  (41)  between  the  measurements  and  trajectory  state, 
we  will  be  required  to  linearize  the  estimation  at  some  stage  of  the 
derivation.  It  is  most  convenient  to  linearize  in  (42)  new  so  as  to 
obtain  the  common  Gauss-Hewton  iteration  sequence.  Let  x°Uk)> 

k=l ,  2,  —  ,n  be  a  guess  trajectory.  We  will  discuss  the  sources  of 
this  guess  later.  Also,  let  6x(tk)  =  x(tk)  -  x°(tk)  and 


H^k) 


3hi (x(tk)) 

3x(tk) 


Then  instead  of  minimizing  (42)  we  will  minimize 


m 


n 


W 


(r.(tk)  -  H.(k)6x(tk)  -  g"|(k)b)  , 


(43) 


where  r.  (tk)  =  z.j(tk)  -  h(x!j(tk)).  VJe  proceed  in  the  usual  way  by 
differentiating  (43)  with  respect  to  5x(tk)  and  b  and  manipulating.  This 
procedure  leads  to  the  set  of  normal  equations, 
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The  matrices  Ak  are  either  3x3  are  6x6,  Ak  n+^  are  3xp  or  6xp 
depending  on  the  dimension  of  the  state  vector  x(tk).  The  vectors 
y(k)  are  either  3-vectors  or  6-vectors.  y(n+l)  is  p-vector  and 
An+i  is  pxp. 

Although  (44)  may  be  of  very  large  dimension  due  to  the  large  number  of 
trajectory  time  points  t^,  the  simple  structure  of  (44)  allows  a  simple 

decomposition  which  leads  to  a  simple  solution.  Let  L  be  a  block  lower 
triangular  matrix  and  let  U  be  a  block  upper  triangular  matrix.  Then  we 
can  find  an  L  and  U  such  that  the  coefficient  matrix  A  in  (44)  can  be 
written  as  A  =  LU.  Equating  corresponding  blocks  of  A  with  those 
of  the  LU  product  gives 


L 
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where 


Pk  =  Ak  Ak,n+1 


Q  =  An+1  *  k?-|Ak,n+lPk 


(52) 


(53) 


The  solution  of  (44)  leads  to  the  solution  of  simple  set  of  equations  which 
reduces  to 

Vx(tk)  =  y(1)  i=l,n  (54) 

Qb  =  y(n+l )  -  l  A^n+1  5x(tk)  (55) 

k-1 


Ak6x(tk)  =  j  H^(k)(r.(k)  -  gI(k)b)/R.(tk)  (56) 


Covariance  estimates  of  5x(tk)  and  b  are  also  easily  computed.  The 
above  linear  least  squares  minimization  is  repeated  with  x!j(tk) 
replaced  with  x?(tk)  +  6x(tk)  until  convergence.  The  batch  processor 
is  incorporated  into  a  BET  as  shown  by  the  following  flow  diagram. 
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PREPROCESSOR 


!  Delete  wild  data  points  and  catastrophically  biased  measurements. 

I  Compute  measurement  noise  variances  and  guess  trajectory. _ 


BATCH  PROCESSOR 

j  Using  all  measurements  compute  improved  trajectory  state  estimate 
i  and  measurement  biases .  _  _ 


+ 


ADAPTIVE,  VARIABLE  LAG  SMOOTHER 

j  Compute  smoothed  position,  velocity,  and 
1  acceleration  states  of  trajectory _ 


+ 


OUTPUT  PROGRAM  j 

Compute  quantities  desired  by  Range  user 
which  are  derivable  from  position,  velocity,  1 
acceleration,  and  air  weather. 


As  shown  in  the  flow  diagram  a  preprocessor  is  used  for  data 
editing  in  the  batch  processor  BET.  The  preprocessor  fits  a  curve 
to  each  measurement  sequence,  deletes  points  which  are  too  far  from  the 
fitted  curve,  and  computes  the  variance  of  the  differences  between  the 
measurements  and  curve  fit  to  use  as  a  processing  variance  in  the  batch 
processor.  The  preprocessor  also  provides  a  guess  trajectory  to  the 
batch  processor  and  interpolates  measurements  to  a  comnon  set  of  times 
if  the  measurements  are  not  synchronized.  The  measurements  need  to  be 
synchronized  in  the  batch  processor  in  order  that  each  type  of  measure¬ 
ment  have  an  influence  on  the  bias  estimate  obtained  for  another  type 
of  measurement.  Thus,  for  example,  if  the  radar  measurements  were 
processed  at  a  different  set  of  times  than  optical  measurements,  the 
optical  measurements  would  not  have  any  influence  on  the  determi nati on 
of  the  radar  bias  estimates  and  vice  versa.  It  is  the  biases  that 
furnish  a  tie  between  times  in  the  batch  processor.  It  is  obvious 
from  the  structure  of  the  batch  processor  equations  given  in  (44)  and 
(46)  that  if  we  have  no  bias  terms  to  estimate  then  the  (n+l)st  row  and 
column  of  blocks  in  (44)  will  be  zero  so  that  the  batch  processor 
reduces  to  an  n(k)-station  solution  at  each  time  point  with  n(k)  is 

the  number  of  instruments  at  t, 

k 

At  each  t.  a  point  x°(t,)  along  the  guess  trajectory  is  computed 
by  the  preprocessor.  If  n^,  nk>_  2,  uptics  stations  are  available  at 

tk’  x°^tk^  is  an  l\'statl’on  optical  solution.  If  n,<  2,  x°(tk)  is 

computed  from  radars  such  that  each  component  of  position  of 

x°(tk)  is  the  median  of  the  correspond!- ng  components  of  the  individual 

radar  cartesian  position.  If  any  measurements  considered  for  the  guess 
are  too  far  from  the  guess  position  solution,  these  measurements  are 
temporarily  removed  from  consideration  and  the  guess  trajectory 
solution  repeated.  If  velocities  are  also  being  estimated  by  the 
batch  processor,  a  guess  solution  for  the  velocity  components  is 
computed  using  the  range  rate  measurements  and  the  cartesian  guess 
positi on. 


Since  the  batch  processor  produces  only  a  trajectory  position 
solution  or  a  position  and  velocity  solution,  an  additional  estimator 
is  required  to  produce  the  derivative  states  for  the  trajectory.  The 
estimator  used  is  an  adaptive,  optimal  variable  lag  smoother  which 
uses  the  raw  batch  processor  state  estimates  to  obtain  smoothed 

estimates  (x$,  xs,‘xs),  (y$,  y$,*ys),  (z$,  z$,  z$).  The  smoothed 

estimates  are  obtained  independently  for  each  coordinate  direction. 
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Let  s ( ) .  a  3-vector,  represent  any  of  three  coordinates  of 
the  trajectory,  i.e.,  s=  (x,  x.’x)  or  s=  (y,  y,'y)  or  s=  (z,  z.’z). 
Assume  that  s ( t^)  obeys  the  discrete  state  equation. 


s(tk+l)  =  *U)s(tk)  +  y(Mu. 


(57) 


It  is  assumed  that  measurements  of  s(tk)  are  equally  spaced  in  time  with 
tk+i*  t^  =  A,  the  measurements  being  a  position  coordinate  or  position  and 
velocity  coordinate  from  the  batch  processor  state  estimate.  Assuming 
constant  acceleration  the  state  transition  matrix  is 


$(A)  = 


1  A 

0  1 

0  0 


A2/ 2 
A 

1 


(58) 


In  (57)  u  is  an  unknown  constant  scalar  forcing  function.  y[t\)  is  the 
vector 


(a)  =  [a  a2/2  a3/6] 


(59) 


In  order  to  avoid  complexity  of  notation,  we  consider  the  smoother 
equations  when  the  batch  processor  is  estimating  only  a  3-vector  of 
position  so  that  a  position  is  the  only  measurement  input  to  the 
smoother.  Let  s(k)  be  the  smoothed  estimate  at  time  tk  given 

measurements  m(t^),  m^), —  »m(tk+n  )•  Th^s,  the  smoother  lags 

by  n^  data  points  where  n^  is  variable.  Let  A(k/k)  be  the  optimal 
estimate  at  time  tk  given  measurements  t^ ,  t^, — .  t^,  i.e.,  the  filtered 
estimate.  Then 

n, 


p}k)  =  P1  (k/k)  +  [  l  ( i  A )  HTH4»  ( i  a  )  ]/r  ( k ) 


(60) 


i  =  l 


s(k)  =  s(k/k)  +  P(k)[  l  *T(iA)HT^n(tk+i )-H$(iA)s(k/k)^]/R(k)  (61) 
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H  =  [1  00],  R(k)  =  measurement  noise  variance 


s(k+l/k)  =  <i>(A)s(k/k)  (62) 

P(k+l/k)  =  $(A)P(k/k)$T(A)  +  q(kb(A)YT(A)  (63) 

P1 (k+l/k+1 )=P  (k+1/k)+HTH/R(k)  (64) 

S (k+l/k+1 )  -  s(k+l/k)  +  P(k+l/k+l)HT(m(k+1)-Hs(k+l/k))/R(k)  (65) 


Equations  (60)  -  (65)  are  not  implemented  directly.  Instead,  these 
equations  are  implemented  in  a  matrix  square,  root  formulation  in 
order  to  ensure  smoother  and  filter  stability.  As  in  the  matrix 
square  root  formulation  of  the  Kalman  filter  BET  all  covariance 
matrices,  P(k),  P(k+l/k),  and  P(k/k)  are  written  as  LU  or  LOL1 
where  L  is  lower  triangular  and  D  is  diagonal.  The  filtering  and 
smoothing  equations  (60)  -  (65)  are  then  expressed  in  terms  of  L 
and  D  and  a  covariance  is  only  computed  in  order  to  output  error 
estimates. 

The  smoother  adapts  to  changing  conditions  by  computing  R(k)  from 
filter  residuals  and  estimating  q(k)  from  the  forward  residual  stack, 

r(k+l)  =  m(tk+^)  -  H$(lA)s(k/k) ,  1=1, n.  The  lag  n^  in  (60)  and  (61) 

is  less  than  the  fixed  integer  n  and  is  a  function  of  q(k);  0<n.<n 

with  equality  n^=  0  for  small  q(k)  and  n^  =  n  for  large  q(k). 


Optimal  Instrumentation  Planning 


The  relative  geometry  of  a  trajectory  and  an  instrumentation  plan 
(IP)  determines  the  quality  of  output  of  a  BET.  It  is  natural  then 
for  someone  who  is  producing  BET's  to  be  concerned  with  the  geometry 
of  the  instrumentation  plan  used  on  a  mission.  It  is  also  a  natural 
extension  of  a  BET  program  to  develop  a  program  to  compute  an  optimal 
instrumentation  plan.  Thus,  as  a  derivative  of  our  batch  processor 
BET  we  have  developed  optimal  instrumentation  planning  programs  for 
DOVAP,  cines,  and  radar.  The  algorithm  for  optimal  instrumentation 
planning  will  be  briefly  described  for  the  radar  case.  The  basic 
algorithm  also  applies  to  cines  and  Dovap  with  modifications  which 
will  be  mentioned. 
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Given  a  set  of  n  time  points  t^  which  entirely  span  a  nominal  flight 
path  specified  by  a  range  user  and  the  corresponding  position  vectors 
x.,  1=1,  n,  to  the  flight  path  a  set  of  M  existing  radar  sites  is 
selected  which  minimizes 

CM  "  f,  “(V  cov(*1> 


The  quantity  may  be  interpreted  as  the  weighted  sum  of  error 

estimates  which  would  result  if  radar  measurements  from  the  M  radars 

were  processed  by  the  batch  processor  BET  program.  The  w  are  a  set  of 

i 

weights  used  to  attach  more  importance  to  some  trajectory  points  than 
to  others.  Cov(x1-)  are  the  3x3  covariance  matrices  of  the  errors  in  the 

estimates  of  the  position  vectors  x^ . 

The  effect  of  changing  the  number  of  radar  sites  IT  to  be  used  in 
instrumenting  the  flight  path  is  easily  computed  using  the  optimal 
site  selection  program.  Using  a  minimum  number  of  sites  M-j  and  a 

maximum  number  os  sites  M2,  the  selection  program  sequentially  obtains 

the  optimal  IP  for  M-j  sites,  i1-|+l  sites,  — ,  M2  sites  so  that  the 

effect  of  additional  sites  on  the  error  estimates  can  be  readily 

determined.  Besides  the  relative  geometry  of  the  radar b  and  flight 
path,  the  selection  program  also  considers  past  radar  performance 
statistics  in  the  form  of  measurement  noise  variances. 


Rather  than  pursuing  a  global  minimum  of  which  would  require 

the  solution  of  a  very  large  combinatorial  programming  problem,  a 
local  minimum  of  is  achieved  through  the  use  of  an  IP  improvement 

algorithm.  The  following  sequence  of  steps  are  executed  in  the  IP 
improvement  algorithm. 


a.  Starting  with  an  initial  IP  having  M  arbitrary  radar  sites 
construct  a  modified  IP  having  M+l  sites  by  adding  the  site  from 
among  ail  remaining  radar  sites  which  results  in  the  smallest 


value  of  C 


M+l  ‘ 


b.  Delete  the  radar  site  from  the  modified  IP  which  results  in 
the  smallest  value  of  Cu. 

c.  Repeat  steps  a  and  b  until  no  further  decrease  in  the  value  of 
in  step  b  can  be  achieved. 


The  algorithm  will  terminate  when  the  radar  added  in  step  a  is  the 
same  radar  deleted  in  step  b.  The  minimum  of  C..  achieved  in  step  b  is 
local  in  the  sense  that  it  is  dependent  on  the  arbitrary  initial  IP 
with  which  the  algorithm  started.  The  existing  sites  from  which  the 
radars  in  the  IP  are  selected  muse  satisfy  some  basic  constraints.  For 
a  radar  to  be  considered  for  selection  the  elevation  angles  cannot  be 
too  small  or  too  large  over  a  large  portion  of  the  trajectory. 

Optimal  instrumentation  planning  programs  have  also  been  wr.tten 
for  DOVAP  and  cinetheodolite.  The  DOVAP  instrumentation  planning 
program  was  the  first  developed  and  was  used  for  all  DOVAP  IP's 
for  some  time  before  DOVAP  was  phased  cut  at  WSMR.  The  use  of  this 
instrumentation  planning  algorithm  resulted  in  an  average  of  25 i 
saving  in  the  number  of  DOVAP  receiver  sites  required  to  achieve  a 
given  quality  of  trajectory  data.  The  application  of  this  algorithm  to 
cine  instrumentation  planning  is  considerably  more  difficult.  The 
difficulty  does  not  occur  in  the  application  of  the  IP  improvement 
algorithm  but  in  implementing  the  several  constraints  which  must  be 
met  before  a  cine  site  can  be  considered  for  inclusion  in  an  IP. 
Constraints  which  must  be  satisfied  for  each  cine  site  considered  are: 
minimum  image  size  readable  on  film,  maximum  tracking  rates,  sun 
angle,  and  flight  safety  evacuation  area.  Development  of  these  optimal 
IP  programs  are  described  in  [22]  -  [25]. 

Directions  of  BET  Development  at  WSMR 

The  development  of  BET  techniques  at  WSMR  is  continual.  As 
mentioned  previously  measurement  data  editing  and  adaptive  filtering 
are  two  difficult  features  of  a  BET  program  which  must  be  successfully 
implemented  to  insure  good  3ET  performance. 

We  have  been  very  active  at  WSMR  in  the  development  of  robust 
estimation  methods  and  their  application  to  several  measurement  editing 
problems  in  data  reduction,  [15J  -  [18].  We  have  been  highly  successful 
in  these  applications  of  robust  estimation  to  measurement  editing  [15]. 

We  have  applied  these  methods  to  the  preprocessor  in  the  MISTE  program, 
to  instrument  calibration  problems,  the  N-station  Davis  solution,  and  are 
currently  developing  some  robust  filtering  methods. 
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Before  describing  the  application  of  robust  estimation  methods  we 
need  to  answer:  What  are  robust  estimation  methods  and  how  do  they 
apply  to  data  editing?  In  answer  to  the  first  part  of  the  question 
we  can  briefly  describe  robust  statistical  methods  as  those  which 
perform  well  under  a  wide  variety  of  underlying  probability  density 
functions  of  the  measurements  or  in  the  presence  of  measuranents  from 
contaminating  distributions.  In  answer  to  the  second  part  of  the 
question, we  are  probably  not  very  concerned  about  the  performance  of 
data  reduction  procedures  under  a  wide  variety  of  underlying  distribution 
functions  for  the  measurements^ but  are  mainly  concerned  about  the  per¬ 
formance  of  our  methods  in  the  presence  of  observations  from  contaminating 
distributions,  i.e.,  outliers  or  wild  data  points.  In  data  reduction  we  are 
usually  interested  in  estimating  the  parameters  in  some  postulated  linear  or 
nonlinear  model  of  a  process.  Thus,  in  data  reduction  we  are  specifically 
interested  in  developing  methods  for  linear  and  nonlinear  regression  which 
are  insensitive  to  a  large  percentage  of  outlying  observations.  The  usual 
methods  of  least  squares,  weighted  least  squares,  maximum  likelihood,  etc, 
used  in  data  reduction  for  estimating  parameters  in  a  model  become  useless 
in  the  presence  of  outliers.  To  quote  Huber  [19],  "even  a  single  grossly 
outlying  observation  may  spoil  the  least  squares  estimate  and  moreover 
outliers  are  much  harder  to  spot  in  the  regression  case  than  in  the 
simple  location  case." 

The  most  popular  and  most  developed  of  the  robust  methods  for 
use  in  parameter  estimation  problems  are  the  M-estimates  of  Huber.  Given 
the  linear  model 
P 

•^i  “  I  x  .  .e .  +  e.  i=l  ,n  (66) 

j=l  J  1 

where  6j  are  unknown  parameters  to  be  estimated.  The  M-estimates  of 
9.j  minimize 

n 

Ip 

i  =  l 


where  p(-)  is  some  suitable  function  and  s  is  a  measure  of  dispersion  of 
the  residuals.  Rather  than  specifying  the  function  r,  M-estimates  are 
usually  specified  in  terms  of  their  *  function  which  is  the  derivative  of 
p.  Several  ip  functions  have  been  proposed  in  the  literature  [15],  but 
we  will  only  consider  a  t \>  function  proposed  by  Hampel  [20].  This  i>  is 
given  by 


'  y .  -  [x..a. 
1  4  ij  j 


(67) 
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a  sgn  (x) 


U'(x)  =  „  x  -  c  sqn(x) 

b  -  c 

0 


a_v;  x ,  <b 


b< I x  <c 


,X;^C 


A  graphical  representation  of  this  function  is  given  below 

<  <L>(X) 


(68) 


X 


.  ~L n  or^er  to  minimize  (67)  we  differentiate  with  resoect  to  e,  where 
0  is  a  p-vector ,  6  =  [o1 ,  e9 — 9n],  This  gives 


where  0  is  the  M-estimate  of  0  and  =  col(x^,  xi2’"'_'  xip^‘ 

(69)  is  the  analog  of  the  normal  equations  in  least  squares 
estimation.  Specifically,  if  <p(x)  is  linear  for  all  x,  we  have  a 
least  square  estimation  process.  (69)  can  be  rewritten  as 


^  wi(0)xl!(yi  -  x.e)  =  0 


where 


w^e)  = 


(npi) 


y.  -  x.f 

i 


(70)  is  solved  iteratively  as  follows.  Starting  at  an  arbitrary  point 
0^  of  the  iteration  sequence,  we  replace  (70)  by 


l  w  (e(k))xj(y.  -  x  e(k+1))  =  0 
i-1  1  iii 


Solving  (72)  for  e^k+1^  gives 


-  l  W.(e^) 
,i=l  J 


-j*))  j, 


(73)  is  just  a  weighted  least  squares  solution.  Thus,  we  apply  a  weighted 
least  squares  algorithm  iteratively  to  obtain  an  M-estimate.  Examining 

(71)  we  see  that  for  large  residuals,  y..  -  x.e>>s,  the  weight  w.  goes  to 

zero.  Thus,  large  residuals  are  weighted  out  of  this  M-estimate  solution. 
The  robust  dispersion  measure  s  which  is  used  is  the  MAD  (Median  of  the 
Absolute  Deviations)  estimate  given  by 

s  =  median|r.|  Z^.6745  (74) 


(74) 


where  =  y.  -  x^a,  Hampel  has  shown  that  the  MAD  estimate  is  the 

most  robust  measure  of  dispersion.  The  above  method  for  compu¬ 
ting  an  M-estimate  is  iterative  and  therefore  requires  a  starting 
solution  to  be  specified.  The  required  closeness  of  the  starting 
solution  to  the  final  solution  depends  on  the  application  and  the  type 
of  ii>  function  used.  Often  an  ordinary  least  squares  solution  provides 
a  sufficiently  good  starting  solution.  In  some  cases  it  is  necessary 
to  use  a  starting  solution  which  is  more  robust,  see  [16]. 

The  application  of  robust  estimation  in  the  BET  preprocessor 
provided  our  original  motivation  for  the  development  and  application 
of  robust  methods  'in  data  reduction.  In  the  preprocessor  we  have 
the  time  history  of  each  measurement  function  for  its  entire  span  of 
observation  on  a  trajectory.  The  preprocessor  divides  this  interval 
of  observation  into  equal  segments  of  T  seconds  except  for  a  final 
segment  either  longer  or  shorter  tnan  T.  Over  each  of  these  segments 
a  polynomial  ,  usually  a  quadratic,  is  fit  to  the  measurements.  Alter¬ 
natively,  a  cubic  spline  might  be  fit  to  the  entire  span  of  measurement 
data  using  the  end  points  of  the  T  second  segments  as  knot  times.  Using 
an  M-estimation  procedure  we  estimate  the  parameters  of  the  polynomial 
and  delete  those  points  which  are  given  zero  weights  in  the  M-estimate 
process  or  whose  residuals  are  greater  than  k*s  where  k  is  a  suitable 
constant.  The  processing  variance  is  computed  as  the  average  squared 
error  of  the  remaining  residuals.  The  following  data  set  is  from  a 
real  preprocessor  application.  The  data  are  a  sequence  of  azimuth 
measurements  from  a  cine. 


Residuals  from 

Observations  Robust  Fit 


Residuals  from 
Least  Squares  Fit 


-1.70987 

.000012 

-.157774 

-1.70942 

-.000004 

-.000204 

-1.70893 

.000003 

.105480 

-1.70845 

-.000015 

.159227 

-1.70793 

-.000010 

.161087 

-1.70741 

-.000021 

.111021 

-1.70682 

.000022 

.009099 

-1.70626 

.000019 

-.144780 

-1.70571 

-.000010 

-.350595 

-1.70510 

.000005 

-.608277 

-1.70449 

.000004 

-.917885 

1.43777 

3.141637 

1.862231 

1.44602 

3.149243 

1.456410 

-1.70257 

-.000007 

-2.158177 

1.44667 

3.146558 

.473139 
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The  residuals  from  the  robust  fit  which  were  obtained  using  a 
Hampel  p  with  breakpoints  (2.5,  5.0,  7.5)  show  exactly  which  points 
are  outliers.  The  least  square  residuals  provide  no  information 
about  the  outl iers. 

Another  application  in  BET  processing  in  which  we  have  successfully 
used  M-estimates  is  in  calibration.  For  example,  a  reasonable  error 
model  for  a  laser  tracker  is 


9,  -  e„tanEc -cosA,. .  -  ectanE„  .sinA„  .  -  e, /cos  E  . 

3  4  sj  sj  b  sj  sg  o  sj 

e?  +  e4sinAsj  -  e5cosA$j 

where  R$ j ,  A$ ^ ,  Eg^  are  the  surveyed  range,  azimuth,  and  elevation  of 

the  jth  target  calibration  target.  We  have  multiple  range,  azimuth,  and 

elevation  observations,  denoted  by  R..,  A..,  E..,  i=l,  N.,  j=l,  m,  for 

U  1J  1J  j 

each  of  the  M  calibration  targets,  then  aR^.  =  R^  -  ,  aA^  =  A.j  -  A^  , 

and  AE.  .  =  E-.  -  E  ..  The  unknown  parameters  e.,  j=l  ,7  are  to  be  estimated. 

'  J  1 J  s  J  J 

The  following  example  illustrates  the  application  of  M-estimates  to  the 
calibration  of  a  laser  tracker.  There  were  approximately  250  observations 
for  each  calibration  target.  The  estimation  used  a  Hampel  ip  function  with 
breakpoints  (2.5,  5.0,  7.5).  The  results  of  this  calibration  are  summarized 
in  the  table  on  the  following  page  by  tabulating  the  number  of  residuals 
for  each  target  lying  in  each  of  the  four  regions  of  the  Hampel  p.  We  only 
show  the  positive  side  of  the  Hampel  ip  with  the  number  of  residuals  in 
each  region  being  the  sum  of  the  numbers  of  residuals  in  the  positive  and 
corresponding  negative  side  of  the  ip  function.  Targets  13-20  are  dumped 
readings  of  the  first  eight  targets.  From  the  table  it  can  be  seen  that 
most  of  the  readings  from  several  target  boards  are  outliers,  particularly 
for  the  dumped  readings.  A  least  squares  estimation  of  the  calibration 
parameters  failed  miserably  on  this  example.  Although  this  example  is 
extreme  for  this  application,  having  about  22%  contamination  by  outliers, 
it  well  illustrates  the  power  of  the  M-estimate  method  in  dealing  with 
outliers. 


AR.  .  = 
ij 

aA..  = 
U 


NUMBER  OF  RESIDUALS 


roOOOOOO©OOO**<XCOOv0OOOO 

vO 


OOOOOOOOOOOOOOOOr-OO 


OOOOOOOcvjr^OOOUDmOO«=rcMvo 
U">  CO  < —  CD  <r  CM 


OCMNOtMMNINJfOTr-NCOtj'tMr-aDOCO 
mLDror-.*cr^rmc\Jc7'COLO. — ■  ro  n  cm 

MNCVjCVICVJC\JCMC\ir-CMMW  CM  CM 


CO*--OOOOOC* 


OO-cfCMCMOvoiOCM^TcO, 
VO  ■VT  ro  CTV 

CM  CM 


OOOOOOOOOr-OLOOOOOOOOO 

ir> 


OOOOOOOOCMOCOiOOOOOOOOO 

«r  m 


O  . —  NOMWNUlNririN^COOJ' —  CO 

mir><nr".«3-«3-oo>'-  ■?  »r  w  n  c  cm 

C\i  CM  CNJ  CM  CM  CM  CVJ  CM  CM  CM  > —  CM  CM  CMCM 


COOOOOOOOS 


OOOOOOOOOOOOOvnOOSOOCTVCO 

CM  CM  r— 


OOOOr—  OOOOIO 


cm  m  vo  O  cr*  vo  o 

O  00  lO  ct> 


lOCMr^O'—  CMr^vncrtavo^-covncMvococMr^t— 

VOtOVOr—CO  CM  CO  ro  CO 

!  CM  CM  CM  CM  CM  CM  CM  CM  CM  CM  * —  r—  . —  f —  r—  r- 


DISTRIBUTION  OF  CALIBRATION  RESIDUALS 


The  robust  M-estimates  of  Huber  can  also  be  developed  to  apply 
to  recursive  filtering.  The  development  of  these  ideas  for  filtering 
are  described  in  [18]  and  [26].  Let  x(t, )  be  the  state  vector  of 
a  linear  dynamic  model 


x(tk+1)  =  <i>(  k)x  ( tk)  +  w(k) 


(75) 


where  t(k)  is  the  state  transition  matrix  and  w(k)  is  an  n-vector  of 
zero  mean  Gaussian  noise  with  covariance  Q(k).  Let  scalar  measurements 

z(tk)  =  H  x(tk)  +  v ( k)  (76) 

be  available  where  v(k)  is  a  random  error  term  which  may  be  contaminated 
by  outliers.  A  robust,  pseudo-maximum  likelihood  filter  for  the  case  is 
given  by  the  equations 


P^i/uHT  ( z(t.x1)  -  Hx(k+l/k+l)\ 


c(k+l/k+l )  =  x(k+l/k)  +  1  t^k+1 

sk+l 


’k+1 


(77) 


Pk+1  =  Pk+l/k  '  •  *' 


/z(tk+1)  -  Hx(k+l/k) 


‘  k+1 


s2  Ip  hthp 
k+1  k+l/k"  MVl/k 


where  <|i‘  is  the  derivative  of  \p  sod  pk+yk  =  $Pk$  +  Q|, 


x(k+l/k)  =  4>x(k/k) 


(79) 


and  we  use  a  robust  measure  of  dispersion  obtained  from  the  residuals 


(77)  is  nonlinear  in  x(k+l/k+l)  and  is  iterated  till  convergence. 

We  are  testing  this  robust  filtering  scheme  on  simulated  and  real 
data  to  determine  its  ability  to  filter  measurement  data  with  large 
amounts  of  contamination  by  outliers.  We  are  also  working  at 
extending  the  above  ideas  to  smoothing. 
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RIDGE  REGRESSION  WITH  UNDERSPECIFIED  MODELS 


1 .0  INTRODUCTION 

w«ll  known  to  analysts  associated  with  regression  studies  are 
the  probleas  related  to  Ill-conditioned  matrices  —  problems 
lending  to  loss  of  precision,  grossly  inaccurate  (inflated) 
estimates  of  the  parameters  and  gross  underestimates  of  the  er¬ 
rors  in  the  estimates.  Inflated  estimates  are  those  which  depart 
unrealistically  (*  3  sigmas)  from  a  priori  estimates.  The 
authors  have  analyzed  four  hundred  experiments  with  real  and 
simulated  data  to  determine  the  circumstances  associated  with 
Inflated  estimates.  The  results  of  the  study  nay  be  summarized 
briefly  and  qualitatively. 

For  purposes  of  discussion  it  is  convenient  to  consider  two 
general  classes  of  problems.  Class  A  is  defined  by  the  presence 
of  a  relatively  small  quantity  of  measurements  accompanied  by 
somewhat  large  random  errors  and  high  correlation  among  the 
errors  in  the  adjusted  parameters.  Class  B  encompasses  a  broad 
family  of  regression  problems  (notably  orbit  determination) 
wherein  the  adjustment  exercise  has  available  an  enormous  re¬ 
dundancy  of  measurements;  the  random  errors  in  the  measurements 
are  relatively  unimportant;  and  the  errors  in  the  estimates  are 
highly  correlated.  Difficulties  with  computational  precision 
were  excluded  from  this  investigation.  Perhaps  the  two  most 
practical  types  of  regression  models  are  correctly  specified  and 
underspecified.  Underspecification  may,  for  example,  result 
from  deficient  knowledge  or  in  some  cases  from  limitations  in  a 
computer  program.  With  correctly  specified  models  and  Class  A 
problems,  some  degree  of  inflation  was  observed  occasionally, 
and  the  occurrence  depended  largely  upon  the  vagaries  of  the 
random  errors  In  the  measurements.  We  were  unable  to  induce  or 
discover  Inflated  estimates  when  correctly  specified  models  were 
Applied  to  Class  B  problems.  Inflation  could  be  Induced  at  will 
in  both  Class  A  nnd  Class  B  problems  with  underspecified  models. 

A  common  situation  resulting  in  rather  dramatic  inflation  is 

charnc ter izod  by  the  triple  correlation;  two  highly  correlated 
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parameters  In  the  specified  model  and  a  third  parameter  (missing 
from  the  model)  whose  errors  correlate  with  those  of  the  other 
two.  In  the  triple  correlation  a  small  error  in  the  unmodeled 
parameter  may  be  magnified  many  times  into  errors  in  the 
ordinary-least-squares  estimates  of  both  the  modeled  parameters. 
Our  premise  regarding  ill-conditioned  normal  matrices  tends  to 
preclude  consideration  of  the  simple  one-to-one  or  double  cor¬ 
relation  between  the  error  in  an  unmodeled  parameter  and  the 
error  in  a  modeled  parameter,  because  this  type  of  correlation 
is  only  accidentally  associated  with  ill-conditioned  matrices; 
however,  it  may  be  noted  tn  passing  that  in  this  circumstance 
the  error  in  the  estimate  is  at  most  comparable  in  size  to  the 
unmodeled  error,  and  if  the  estimate  is  "inflated",  then  the  un- 
raodeled  error  1b  so  large  as  to  be  generally  detected  and  removed 
or  else  modeled.  In  any  case,  the  simple  one-to-one  correlation 
between  a  modeled  error  and  an  unmodeled  error  results  in  an 
ordinary-least-squares  estimate  which  is  relatively  unchanged  by 
subsequent  ridge  regression. 

One  of  the  most  fruitful  approaches  to  the  problem  of  inflated 
estimates  was  developed  by  Hoerl  and  Kennard  ( 6  j  [  7]  and  labeled 
by  them  "ridge  regression."  Hoerl  and  Kennard  point  out  an  im¬ 
portant  deficiency  in  ordinary-least-squares,  point-estimation 
procedure:  In  the  case  of  high  correlation  among  the  errors  in 

the  parameter  estimates,  there  may  be  a  gross  inflation  of  the 
adjustment  vector  in  ordor  to  achieve  a  final  minuscular  reduc¬ 
tion  in  the  sum  of  squares  of  the  residuals.  The  Hoerl-Kennard 
((IK)  estimation  process  is  inherently  Bayesian  in  nature.  It 
assumes  expected  values  of  zero  for  the  adjustable  parameters 
and  tends  to  constrain  the  adjusted  values  as  close  as  possible 
to  zero  without  unduly  enlarging  the  measurement  residuals.  The 
HK  estimator  is  biased  but  offers  potential  reduction  in  mean 
square  error.  The  mathematics  supporting  the  HK  estimator  is 
bnRed  upon  a  known  fixed  regression  model  and  classical  un¬ 
weighted  least  squares.  The  HK  estimator  is  not  compatible  with 
standard  Bayesian  minimum  variance  regression  programs  such  as 
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those  associated  with  orbit  detormi  nation .  In  this  paper  we 
develop  and  discuss  the  properties  of  an  unbiased  ridge  estimator 
with  application  to  problems  where  the  regression  model  is  under- 
specified.  The  resulting  estimator  is  readily  compatible  with 
standard  Bayesian  minimum  variance  regression  computer  programs. 

For  purposes  of  comparison  it  will  be  helpful  to  discuss  the 
standard  estimator. 

2.0  STANDARD  BAYES1 AM  MINIMUM  VARIANCE  ESTIMATOR 

Wr  assume  the  standard  linear  approximation  to  the  general  non¬ 
linear  multiple  regression  problem: 

Y  =>  X  P  -f  €  (2.1) 

111 

where  Y  is  (nxl)  and  denotes  the  measurement  vector;  X  is  (nxj) 
and  of  rank  J  and  represents  the  partial  derivative  matrix  of 
nonstochastic  elements  relating  the  mean  values  of  the  measure¬ 
ments  to  the  adjusted  parameters;  0  is  (jxl)  and  designates  the 

true  fixed  but  unknown  parameter  vector;  €  is  (nxl)  and  con- 

m 

sti Lutes  the  measurement  error  vector.  We  assume  E(e  )  is  zero 

m 

and  K(e  c  ')  =  VAR(e  )  =  E  ,  where  E  is  (nxn)  and  known, 
m  a  mm  ro 

In  addition,  we  have  prior  information  consisting  of  a  k  element 
parameter  vector  0Q ,  which  unblasedly  estimates  R0  and  a  K  element 
parameter  error  vector  c  .  is  known  from  introspection  or  from 

previous  independent  measurements.  Therefore 

Po  =  R0  +  ep  (2.2) 

where  K,  originally  suggested  by  Thiel  [16],  is  (kxj)  of  rank  k 

and  consists  of  known  nonstochastic  elements.  If,  for  example, 

R  •  (l  0],  where  I  is  a  (kxk)  unit  matrix  and  0  is  a 

[k  x  (i-k)]  zero  matrix,  then  0Q  represents  estimates  of  the 

Mrsi  k  elements  of  (3.  Equation  (2.2)  assumes  0Q  is  random  and 

hence  represents  a  departure  from  the  Bayesian  approach,  which 

assumes  a  prior  distribution  on  0,  here  considered  fixed.  In 

addition,  E(«‘  )  is  zero  and  E(e  e  ')  =  VAR(e  )  =  E  ,  where  E  is 
P  PP  p  p  P 


(kxk)  and  known.  Furthermore,  we  assume  COV  (fm>  fcp)  zero. 

In  order  to  include  the  prior  information  in  the  estimation  of  0, 
we  combine  equations  (2.1)  and  (2.2)  as  follows: 


-  — 

—  - 

—  — 

Y 

X 

e 

ra 

= 

0  + 

0 

R 

e 

o 

1 

P  _ 

(2.3) 


or  in  an  obvious  change  of  notation 

Y  =  X  £3  +  6 

where 


m 


VAR  ( O 


0 

E 


=  E 


(2.4) 


The  prior  information  has  thus  assumed  the  role  of  measurements. 
Applying  generalized  least  squares  to  (2.4),  we  obtain  the 

A 

following  relation  for  the  estimator  0  of  the  parameter  vector  p: 

0  -  (X'  E_1  X)~A  (X'  E-1  Y)  .  (2.5) 

This  converts  by  simple  substitution  to 

0  (X'  E  ~a  X  +  R'  Z'1  R)"1  (X  '  E  ”1  Y  +  R  '  E  ~ 1  0  )  (2.6) 

m  p  m  p  o 

or  in  an  obvious  change  in  notation 

5  -  P  (X'  Em_1  Y  +  R'  Ep_1  Pc)  .  (2.7) 

Application  of  the  expectation  operator  to  Equation  (2.6)  gives 
K(0)  0,  and  hence  0  is  an  unbiased  estimator.  Application  of 

the  law  of  covariance  propagation  gives 


?A8 


(2.8) 


VAR(0)  =  P  . 

For  computational  convenience,  Equation  (2.7)  may  be  linearized 
to  obtain  the  following  iterative  form: 

A0  -  P  TX'  Em”1  (AY)  +  R'  Ep"1  (A0o>j  ,  (2.9) 

A 

where  A0  is  the  vector  of  corrections  to  the  current  estimates 
of  the  adjusted  parameters;  AY  is  the  vector  of  measurement 
residuals;  and  aP-  is  the  vector  of  differences  between  current 
and  a  priori  estimates  of  the  parameters.  It  can  be  shown  that 

A 

3  is  the  best  linear  unbiased  estimate  in  the  sense  that  it  has 
the  smallest  variance.  Also,  if  the  measurement  errors  are 

A 

normally  distributed,  then  0  has  minimum  variance  among  all  un¬ 
biased  estimates.  Equation  (2.9)  is  widely  used  in  orbit 
determination  programs  and  elsewhere. 

3.0  RIDGE  ESTIMATOR 

Highly  correlated  errors  in  the  parameter  estimates  result  in 

poor  conditioning  of  the  (X  '  E  X  +  R '  E  R)  matrix  in 

m  p 

Equation  (2.0).  The  poorer  the  condition  of  this  matrix,  the 

A 

greater  the  expected  discrepancy  between  0  and  the  true  vector 
0.  On  the  other  hand,  the  worse  the  conditioning,  the  less  is 
the  sensitivity  of  the  residual  sum  of  squares  to  small  de- 

A 

partures  from  0.  Following  the  concept  of  Hoerl  and  Kennard,  we 
impose  an  accessory  condition  upon  the  least-Bquares  criterion, 

A 

thereby  restraining  the  vector  (0  -  0  )  without  greatly  in¬ 
fluencing  the  residual  sum  of  squares. 

Let  R  be  any  estimate  of  the  vector  0.  Then  the  sura  of  squares 
of  the  weighted  measurement  residuals  is  given  by  (Y  -  XB)  *  E 
* (Y  -  XB) .  Using  Lagrangian  constraints,  we  minimize 

¥  -  (B-£>0)  Ep_1  (B-0o)  +  (1/h)  J” (Y  -  XB)  '  E-1  (Y  -  X8)  -  *j  (3.1) 

is  the  multiplier  and  ♦  is  the  total  sum  of  squares. 


where  (1/h) 


We  obtain,  therefore, 


-  —  U  »  E  _1  (B-C  )  -  (1/h)  X'  L-1  (Y  -  XU)  . 

An  p  ° 


(3.2) 


This  reduces  to 

A 

b  -  a*  = 


X  '  E"1  X  +  (h+l)  R'  E  “1  R  1 
m  P 


-1 


•  ‘  X  '  E  -)  Y  +  (h+l)  R'  E  ”1  fl  1  ;  h  >  0  (3.3) 

L  ra  p  oj 


or  in  an  obvious  change  of  notation 


B*  *  Q  X'  E  _1  Y  +  (h+l)  R'  E  “1  9 
m  p  o 


h  >  0 


(3.4) 


Application  of  the  expectation  operator  to  Equation  (3.3)  gives 

A  A 

V.  (|  *)  -  p,  and  hence  0*  is  an  unbiased  estimator.  Application 
of  the  law  of  covariance  propagation  gives 


VAR  (B* ) 


Qf  X  '  E  ~1  X  + 

L  ra 


(h+l)2  R' 


A  A  A  A 

Clearly,  when  h  -  0,  ft*  =  3  and  VAR(3*)  =  VAR(B)  . 


(3.5) 


A  A 

B*  is  related  to  3  as  follows: 

3*  =  QP-1  5  +  Q  h  R'  Ep_1  pQ  .  (3.6) 

linearizing  Equation  (3.3)  results  in  the  following  iterative 
!  orm ; 


Qi  x 


-l 


m 


(AY)  +  (h+l)  R'  E 


-1 


(A9o) 


(3.7) 


where  /<0*  is  the  vector  of  corrections  to  the  current  estimates 
of  the  adjusted  parameters.  Equation  (3.7)  is  the  same  as 
Equation  (:’.'»)  with  the  exception  that  the  input  constant 
(h+ ! )  E(  1  j  has  replaced  the  Input  constant  (Ep-1)  and  hence 

Equal i  n  (3.7)  is  compatible  with  a  broad  class  of  regression 
prog  rams . 
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Subtracting  Equation  (2.8)  from  Equation  (3.5)  gives,  not  un¬ 
expectedly,  a  positive  definite  matrix  for  h  >  0,  and  therefore 

A  A 

0*  with  h  >  0  is  less  "efficient  ’  than  0  for  application  with 
exact  regression  models.  Consider,  however,  application  of  the 

A 

0*  estimator  to  underspecified  regression  models. 

Suppose  we  fit  the  model  in  Equation  (2.1)  to  data  which  obey 
the  model 


=  X  a  +  Z  y  +  €_ 


(3.8) 


whore  Z  is  (nxs)  and  of  rank  s  and  represents  the  partial  deriva¬ 
tive  matrix  of  nonstochastic  elements  relating  the  mean  values  of 
the  measurements  to  the  unadjusted  parameters;  y  is  (sxl)  and 
designates  a  true  fixed  but  unknown  parameter  vector.  The  vector 
8  will  be  estimated  using  the  estimator  given  in  Equation  (3.3). 
Applying  the  expectation  operator,  we  now  obtain 


E(0«)  -  0  +  lx'  E~l  X  +  (h+1)  It'  E~l  r]  X  E'1  Z  y  , 

ra  p  j  m 


(3.9) 


and  hence  0*  is  biased.  VAR(0*)  is  still  given  by  Equation  (3.5), 
and  the  total  mean  square  error  —  sum  of  variances  and  squared 

biases  —  is  given  by 


MSE ( 0* ) 


Trncc  VAR(0*)  +  ^E(0*)  -  0 j '  j^E(0*)  -  0j 


(3.10) 


which  for  h  0,  is  identical  with  MSE((1)  for  the  same  under- 

A 

specified  model.  In  investigating  the  behavior  of  MSE(0*),  we 

/A 

need  the  partial  derivative  of  MSE(0*)  with  respect  to  h. 
because  the  partial  derivative  of  Trace  VAR(8)  is  the  same  as 

A 

the  trace  of  the  partial  derivative  of  VAR(0*),  we  take  the 
partial  derivative  of  VAR(0*)  in  a  straightforward  manner  and 
note  that  it  is  a  continuous  function  of  h  and  is  symmetric 
posit  i vo  definite  for  all  h  >  0.  (It  is  a  null  matrix  for  h  =  0.) 

A 

Consequently,  the  trace  of  VAR(ft*)  is  monotone  Increasing  with 
»< loaning  h  and  in  fnct  approaches  infinity  as  h  approaches 
infinity.  The  derivative  of  j^Fft1*)  -  It  '  ^E(|t*)  -  fjj  is  also 


continuous,  and  tho  matrix  of  the  quadratic  form  of  the  deriva¬ 
tive  is  symmetric  negative  definite  for  all  h  >  0,  approaching 

the  null  matrix  as  h  approaches  infinity.  Therefore 

r  I  ,  r  ~  1 

^  E  ( 0  * )  -  0  I E(0*)  -  0  |  is  monotone  decreasing  with  increasing 
h  but  can  never  be  negative.  Hence,  we  have  a  situation  where 

A  A 

n  t  h  0,  MSE(P*)  -  MSE(fi),  and  at  very  large  values  of  h, 

/V  A 

MSE ( 0* )  >  MSE ( 0) .  We  inquire  whether  there  is  any  intermediate 

A  A 

positive  value  of  h  which  will  result  in  MSE(P*)  <  MSE(8). 

A 

Apparently  there  is,  because  the  partial  derivative  of  VAR(1*) 
evaluated  at  h  =  0  is  the  null  matrix,  whereas  the  partial  deriva- 


il' 


r-  ✓n 

E(P*)  -  f)j  evaluated  at  h  =  0  is  negati 


ve 


live  of  E(P*)  -  ft  | 

L  J  _ 

definite.  Thus  for  some  0  <  h  <  «,  iviSE(0*)  is  a  minimum  and  less 

A 

than  MSEC).  The  authors  have  not  attempted  to  determine  ex¬ 
plicitly  the  value  of  h  corresponding  to  this  minimum  because  the 
value  will  depend  upon  y  which  is  presumably  completely  unknown 
and  unknowable.  Rather  we  adopt  the  empirical  Ridge-Trace 
npproach  used  by  Hoerl  and  Kennard;  i.e.,h  is  chosen  to  have 
that  value  at  which  the  parameter  estimates  appear  to  have 
reached  stability. 


1 . 0  APPLICATION  IN  ORBIT  DETERMINATION 

The  interested  reader  will  find  a  brief  discussion  of  orbit 
determination  in  Appendix  A.  In  order  to  better  illustrate  the 
performance  of  the  ridge  est  mator  with  an  underspecified  model, 
wo  have  chosen  to  present  a  simulation  exercise  rather  than  a 
treatment  of  real  data.  This  application  involves  a  standard 
satellite  orbit  determination  (Cowell,  special  perturbations, 
batch  processing)  in  which  the  adjustable  parameters  include 
orbital  elements  (initial  conditions)  and  radar  coefficients. 

The  orbital  elements  define  the  position  and  velocity  of  the 
satellite  at  epoch.  The  measurements  are  radar  track  data: 
range,  azimuth  and  elevation. 

The  radar  track  data  are  characterized  by  certain  errors  which 
nay  be  expressed  as  linear  terms  in  the  so-called  radar  measure¬ 
ment  equations.  The  radar  measurement  equations,  truncated  so 


as  to  contain  only  terms  of  present  Interest,  are  as  follows: 


A 

z 

*t  +  »! 

+ 

sec  E^ 

measurement 

true  zero  set 

colliraation 

+ 

u  sin  At  tan  Et  -  v 

cos 

At  tan  Et 

mislevel 

+ 

a3  tan  E^ 

+ 

'  A 

non-orthogonality 

random  error 

E 

3 

Et  +  °i 

+ 

e2  cos  Et 

measurement 

true  zero  set 

droop 

+ 

u  cos  At  +  v  sin  At 

+  e 

3  Ctn  Et 

nlslevel  residual 

refraction 


+  e„  (4.2) 

E 

random  error 

Equations  (4.1)  and  (4.2)  plus  the  equations  of  motion  of  the 
satellite  constitute  the  exact  model.  The  underspecified  model 
does  not  Include  terms  in  non-orthogonality  and  residual  re- 
iraotion  In  Equations  (4.1)  and  (4.2),  A  represents  azimuth; 

E,  elevation.  The  zero-set  errors  are  constant  bias  or  offset 
values.  Colllmntion  represents  the  lack  of  perpendicularity 
between  the  radar  beam  and  the  elevation  axis.  Non-orthogonality 
denotes  the  lack  of  perpendicularity  between  the  azimuth  and 
.  ievation  axes.  Mislevel  represents  the  tilt  of  the  azimuth 
plane  —  u  being  the  northward  component  and  v  being  the  east¬ 
ward  component.  This  tilt  is  defined  with  respect  to  the  local 
horizontal  to  the  geodetic  spheroid.  Droop  represents  the  sag 
of  the  radar  beam  axis.  Residual  refraction  refers  to  the  error 
remaining  after  approximate  corrections  have  been  made  for  the 


bending  of  the  radar  beam  traversing  the  atmosphere.  The  random 
errors  represent  noise  in  the  data  and  have  zero  means. 


The  primary  purpose  of  this  experiment  is  to  estimate  the  radar 
coefficients.  In  order  to  accomplish  this,  it  is  necessary 
incidentally  to  determine  the  orbit.  Determining  the  orbit 
constitutes  estimating  the  orbital  elements  at  epoch.  Since 
the  orbital  elements  correlate  only  weakly  with  the  radar  co¬ 
efficients  in  the  geometrical  environment  we  have  chosen,  we 
shall  not  apply  the  principles  of  ridge  analysis  specifically 
to  the  elements,  although  the  elements  will  be  estimated  each 
time  the  radar  coefficients  are  estimated.  The  true  values  of 
the  radar  coefficients,  the  a  priori  estimates  of  these  values 
and  the  uncertainties  (standard  deviations)  in  the  prior  esti¬ 
mates  are  given  in  Table  I. 


TABLE  I 

NUMERICAL  VALUES  OF  RADAR  COEFFICIENTS 


RADAR 

COEFFICIENTS 

TRUE  VALUE  (S  ) 

X 

A  PRIORI 
ESTIMATE  (3qjl) 

S.D.  (aoi) 

*1 

(deg) 

+0.003 

0 

0.003 

a2 

(deg) 

-0.003 

0 

0.003 

u 

(deg) 

-0.002 

0 

0.002 

V 

(deg) 

-0.003 

0 

0.002 

(deg) 

+0.005 

0 

0.003 

(deg) 

+0.002 

0 

0.002 

"j 

(deg) 

-0.002 

— 

— 

°3 

(deg) 

+0.002 

— 

— 
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The  mathematical  adjustment  procedure  1b  iterative  and  is  based 
upon  liquation  (3./),  the  equations  of  motion  ol  the  satellite, 
and  the  radar  measurement  equations.  Initially,  with  h  =  0, 
the  procedure  is  the  standard  one  in  orbit  determination.  After 
a  converged  solution  with  h  =  0  has  been  obtained,  then  ten  or 
so  additional  solutions  with  increasing  values  of  h  are  computed 
in  order  to  define  the  curves  comprising  the  ridge  trace.  A 
pertinent  sub-set  of  the  correlation  matrix  with  h  =  0  is  shown 
in  Table  II  with  elements  rounded  to  three  digits. 

TABLE  II 

CORRELATION  COEFFICIENTS 


n  I  -;o  show  the  root-mean-square  of  the  weighted  measurement 
res  i  dun  Is  (t>*)  vs  h.  The  symbol  <j  is  used  for  the  a  priori 
standard  deviation  in  B  where  the  subscript  1  designates 
1  Ik  i**1  element  in  the  parameter  vector.  Figure  1  shows 
graphically  the  discrepancies  between  the  Individual  estimates 
(0t)  and  the  n  priori  individual  estimates  (ft  i).  Since  this 
i  a  simulation  exercise,  we  may  also  present  graphically  the 
discrepancies  between  the  estimates  (Bt)  and  the  true  values  of 
tin  parameters  In  Figure  2,  therefore,  we  plot 

o,  1 

~  V8  h- 
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I t  Is  observed  that  the  weighted  rms  in  the  measurement 
residuals  (<>*)  has  increased  less  than  1%  as  h  has  increased 
l rom  0  to  30,  the  point  at  which  stability  In  the  estimates  is 
effectively  reached.  The  Inflation  characteristic  of  the 
original  estimates  of  a^  a2,  e^  and  e2  has  been  largely  re¬ 
moved  at  h  =30.  Since  the  errors  in  the  estimates  of  u  and  v 
are  relatively  uncorrelated  with  each  other  or  with  those  of 
any  other  parameter,  their  estimates  are  relatively  unaffected 
by  changes  in  the  value  of  h.  Nof*  so  obvious  in  the  curves  is 
the  fact  that  the  estimate  of  e  corresponding  to  standard 
least  squares  (h=G)  not  only  was  grossly  in  error  but  also  h&u 
the  wrong  sign. 

Composites  of  Figures  1  and  2  are  shown  in  Figure  3,  in  which 

A  A 

plots  of  root-mean-square  of  (ft?  -  ft  ,)/o  .  and  (B?  -  ft.  )/cr 

5  Ol  of  1  1  Ol 

are  given  as  a  function  of  h.  It  is  particularly  encouraging 
that  not  only  do  these  curves  show  marked  reduction  in  dis¬ 
crepancies  as  a  result  of  ridge  regression  but  also  the 
discrepancies  relative  to  true  values  are  smaller  than  dis¬ 
crepancies  relative  to  a  priori  estimates. 

5 . 0  CONCLUDING  REMARKS 

At  first  glance  the  reader  might  be  alarmed  at  the  rather  large 
value  of  30  arrived  at  for  h  in  the  numerical  example.  It 
appears  that  the  prior  information  has  been  given  (nearly)  full 
weight.  Actually  this  is  noi  the  case.  If  the  prior  informa¬ 
tion  had  been  given  full  weight,  then  the  curves  In  Figure  1 
would  show  a  general  tendency  to  be  tangent  to  the  zero  line 
at  h  =  30,  whereas  most  of  them  show  a  strong  disinclination 
to  approach  zero  even  at  h  =  100.  Furthermore,  if  h  had  been 
increased  to  the  point  where  prior  Information  was  given  full 
weight,  the  residual  sum  of  squares  would  —  except  in  a 
prohibitively  unlikely  coincidence  —  have  shown  a  marked  in¬ 
crease.  It  should  also  be  recalled  that  the  prior  information 
entered  the  adjustment  process  in  the  role  of  measurements 
(only  six),  and  in  spite  of  the  large  weight  attached  to  them 
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at  h  .10,  they  could  not  dominate  a  solution  which  included 
over  three  thousand  other  measurements.  In  reality,  the  ridge 
estimation  procedure  has  a  significant  effect  only  upon  the 
parameter  estimates  whose  errors  ate  mutually  correlated,  and 
with  these  the  adjustments  are  minimized  and  portioned  out  in¬ 
versely  according  to  their  a  priori  variances  so  far  as  possible 
without  unduly  enlarging  the  residuals. 
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APPENDIX  A 


NOTES  ON  ORBIT  DETERMINATION 

Preliminary  orbit  determination  uses  only  a  few  observations  and 
is  usually  based  upon  simple  geometrical  and  dynamical  principles 
of  Kepler's  laws.  Historically  this  type  was  developed  to  keep 
track  of  various  planets  and  comets  sufficiently  well  to  point 
telescopes  and  to  recognize  the  objects  on  their  reappearance. 
Preliminary  orbit  determination  in  some  form  is  also  usually 
necessary  to  obtain  approximate  initial  conditions  for  starting 
the  final  orbit  determination.  Escobal  [3]  gives  a  complete  and 
modern  treatment  of  the  mathematics  of  the  various  methods  of 
preliminary  orbit  determination.  Huseonlca  [8]  evaluates  the 
methods  of  [3]  in  terras  of  limits  of  applicability,  accuracy, 
storage  requirements  and  computation  time. 

Final  orbit  determination  makes  use  of  considerable  redundant 
data  and  is  basically  nothing  more  than  well-known  regression 
analysis.  See,  for  example,  Solloway  [13].  In  the  so-called 
batch"  process  all  the  data  are  fitted  at  once,  and  a  single  set 
of  initial  conditions  is  estimated  along  with  a  few  dozen  other 
parameters  relating  to  the  environment  and  the  tracker  charac¬ 
teristics.  Equation  (2.9)  is  characteristic  of  the  differential 
correction  process.  The  steps  are  briefly  as  follows: 

1.  An  orbit  is  generated  using  the  approximate 
initial  conditions  and  the  equations  of  motion. 

2.  Tho  requirod  partial  derivatives  are  obtained 
by  various  means  —  typically  by  analytical 
methods,  variational  equations,  or  finite 

di f lerences. 

3 .  At  time  [joints  coinciding  with  actual  observa¬ 
tions,  position  points  along  the  generated 
trajectory  are  transformed  to  computed 
observations . 
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1.  l)i  1’ lerencing  the  actual  and  computed  observations 
provides  residuals. 

5.  Information  from  Steps  2  and  4  when  supplied  to 
Equation  (2.9)  permits  an  improved  estimate  of 
the  initial  conditions  and  the  other  parameters, 
and  the  process  starts  over  at  Step  1. 

6.  These  iterations  continue  until  some  convergence 
criterion  is  satisfied.  Modern  computer  programs 
contain  a  flexible  bounding  procedure  to  assist 
in  achieving  a  converged  solution. 

The  final  orbit  determination  makes  use  of  a  wide  variety  of  orbit 
generators.  If  numerical  integration  of  the  acceleration  equa¬ 
tions  is  used,  the  process  is  called  "special  perturbations  ' . 

If  analytical  integration  of  series  expansions  of  the  perturbative 
accelerations  is  used,  the  process  is  "general  perturbations". 

Cowell's  special  perturbation  technique  involves  the  direct, 
step-by-step  integration  of  the  total  acceleration,  central  as 
well  as  perturbative,  of  the  satellite.  Encke's  method  differs 
from  Cowell’s  in  that  differential  accelerations  or  the  deviations 
from  a  two-body  reference  orbit  rather  than  the  total  accelera¬ 
tions  are  integrated.  The  var iatlon-of-parameters  method  differs 
from  Encke's  in  that  there  is  a  continuous  rectification  of  the 
reference  orbit.  All  three  of  these  special  perturbation  methods 
are  described  lucidly  by  Baker  [1].  At  the  present  time,  using 
i  li<'  most  refined  techniques  in  special  perturbations,  earth 
•nlolllte  orbits  are  frequently  computed  to  an  rms  position  error 
<>f  only  thirty  feet  relative  to  the  center  of  the  earth  in  a 
single  fitting  of  up  to  sixteen  days  of  tracking  data. 

General  perturbation  schemes  are  particularly  valuable  for  pre¬ 
dicting  over  long  time  periods.  Their  advantages  are  reduced 
computer  time  and  the  fact  that  they  allow  a  clearer  interpreta¬ 
tion  of  the  sources  of  the  perturbative  forces.  Some  of  the 


netter  known  general  perturbations  schemes  are  Brouwer's  [2), 
Kozai'B  [10],  Musen's  [11 ]  and  Frazer’s  [5],  There  are  many 
others.  The  accuracy  of  general  perturbation  methods  suffers  to 
some  extent  because  these  methods  are  all  limited  to  a  very 
simple  description  of  the  earth's  gravity  field.  They  are,  how¬ 
ever,  very  widely  used  because  of  compensating  advantages. 

One  of  the  more  recent  developments  in  orbit  determination  is 
sequential  processing,  typified  by  the  use  of  the  Kalman  filter 
[9].  This  procedure  gives  an  orbit  update  as  each  new  measure¬ 
ment  is  received.  In  accuracy  it  is  essentially  equivalent  to 
the  batch  process.  It  has  been  used  quite  successfully  in  con¬ 
cert  with  the  Encke  process  [14],  In  computational  efficiency 
it  is  superior  to  the  batch  processor  if  frequent  intermediate 
answers  are  required,  but  it  is  less  efficient  if  only  a  single 
final  answer  is  required. 

Finally,  there  are  the  sequential-batch  processes,  which  take 
many  forms  and  applications.  Originally  proposed  by  Swerling 
[15],  the  sequential-batch  process  can  give  normal  weight  to  all 
prior  information,  or  it  can,  by  reducing  the  weight  on  prior 
information,  act  like  a  fading  memory  filter.  Small  randomly 
varying  parameters  may  be  excluded  from  the  estimation  vector  or 
may  be  assumed  to  be  piece-wise  constant  and  estimated  [4], 

Characteristic  of  all  Of  these  orbit  determination  methods  is  an 
extraordinary  number  and  variety  of  transformations.  Reference 
[12]  provides  a  treatment  of  many  of  these  transformations  along 
with  details  of  orbit  generators,  coordinate  systems,  equations 
of  motion,  special  and  general  perturbations,  earth  gravitational 
field,  variational  equations,  sequential  processing,  error 
analysis,  geodetic  datums,  systems  of  time,  etc. 
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